Setup

knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
# Set options for analysis
shutdown = FALSE

options(dplyr.print_max = 100, 
        brms.backend = 'cmdstan',
        brms.file_refit = 'on_change' # "on_change" or "never"
        )
# functions for kable printing. 
print_df <- function(df, width = "auto") {
  kable(df) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "responsive"),
      full_width = F,
      position = "left",
      fixed_thead = T
    ) %>%
    column_spec(2:ncol(df), width = width) %>%
    row_spec(1:nrow(df), extra_css = "white-space: nowrap;") %>%
    scroll_box(width = '100%', height = '500px')
}

packing_ind <- function(df) {
  pack_rows(df,
    "Within-Person Effects", 2, 6) %>%
  pack_rows(
    "Between-Person Effects", 7, 10) %>%
  pack_rows(
    "Random Effects", 11, 13) %>%
  pack_rows("Additional Parameters", 14, 18)
}

packing_apim <- function(df) {
  pack_rows(df,
    "Within-Person Effects", 3, 12) %>%
  pack_rows(
    "Between-Person Effects", 13, 20) %>%
  pack_rows(
    "Random Effects", 21, 26) %>%
  pack_rows("Additional Parameters", 27, 31)
}


# Function to prepare for export
prepare_for_export <- function(model) {
  rep <- model
  rep$coefs <- rownames(rep)
  final <- cbind(rep$coefs, model)
  return(final)
}

Data Preparation

Scales and Scores

predictors_to_center <- c("persuasion",
                          'pressure',
                          'pushing', 
                          'weartime',
                          'barriers',
                          'support',
                          'comf',
                          'reas')


predictors_not_to_center <- c('day', 'userID', 'coupleID', # independent variables that don't need to be centered.
                              'aff', # outcomes (not centered)
                              'pa_sub',
                              'pa_obj', 
                              'anticipAff', 
                              'reactance', 
                              'plan',
                              'isWeekend',
                              'studyGroup',
                              'got_JITAI',
                              'skilled_support'
                              )



all_variables <- c(predictors_not_to_center, predictors_to_center)


df <- df %>% 
  mutate(ss_affect1 = ss_affect1 + 1, # re-code items from 0-5 to 1-6 scale.
         ss_affect2 = ss_affect2 + 1, 
         ss_affect3 = ss_affect3 + 1, 
         ss_affect4 = ss_affect4 + 1, 
         
         aff = (ss_affect1 + ss_affect3 + (7 - ss_affect2) + (7 - ss_affect4)) / 4, #Affect Scale

         
         pa_sub = ifelse(ss_pa == 0, 0, ss_pa_min_solo + ss_pa_min_collab),
         
         pa_obj = ifelse(wear_minutes > 600, minutes_mvpa_non_filtered, NA),
         
         day = day/54,
         
         barriers = (ss_barr_1 + ss_barr_2 + ss_barr_3 + ss_barr_4 + ss_barr_5 + ss_barr_6 + ss_barr_7) / 7,
      
         plan = ifelse(
           is.na(ss_pa_no), 
           ifelse(ss_pa_yes_0 == 1 | ss_pa_yes_1 == 1 | ss_pa_yes_2 == 1 | ss_pa_yes_3 == 1, 1,0),
           ifelse(ss_pa_no ==1, 1, 0)
         ), 
         
         got_JITAI = ifelse(((userID %% 2 != 0) & (trig_target_slot_0 %in% c("A", "C") | trig_target_slot_1 %in% c("A", "C") | trig_target_slot_2 %in% c("A", "C"))) |
                              ((userID %% 2 == 0) & (trig_target_slot_0 %in% c("B", "C") | trig_target_slot_1 %in% c("B", "C") | trig_target_slot_2 %in% c("B", "C"))), 1, 0),
         
         got_JITAI = ifelse(is.na(got_JITAI), 0, got_JITAI),
         
         skilled_support = ifelse(
           (studyGroup == 3 & ((day * 54 + 1) > 28)) | (studyGroup != 3 & ((day * 54 + 1) > 7)), 
           1, 0
           )
         
         )%>% 
  rename(persuasion = sp_psc_more,
         pressure = sp_nsc_more,
         pushing = sp_push_plan,
         anticipAff = ss_affect_con,
         reactance = ss_reactance,
         weartime = wear_minutes,
         # aff_reactance_persuasion = sp_persuasion_aff,
         #aff_reactance_pressure = sp_nsc_aff, # explorative
         #aff_reactance_pushing = sp_push_plan_aff, # explorative
         
         support = sp_emo_pleasure,
         comf = sp_emo_comf,
         reas = sp_emo_reass
         ) 
# pushing can only occur if there was a plan! Therefore, we set it to zero if there was no plan. 
df$pushing[df$plan == 0] <- 0

center Data

prepare_df <- function(df) {
  df_wide <- df %>%
    # center
    bmlm::isolate(by = "userID",
                  value = predictors_to_center,
                  which = "both") %>%
    # distinguish
    group_by(userID) %>%
    mutate(person_mean_pa_obj = mean(pa_obj, na.rm = TRUE)) %>%
    ungroup() %>%
    
    group_by(coupleID) %>%
    mutate(
      role = case_when(
        person_mean_pa_obj == max(person_mean_pa_obj, na.rm = TRUE) ~ 'A', 
        TRUE ~ 'B'
    )) %>%
    ungroup() %>%
  
    # Reshape the data to wide format
    pivot_wider(
      id_cols = c(coupleID, day),
      names_from = c(role),
      values_from = -c(1, 2, 4)
    )
    
    
  # Double Stacking Outcomes 
    
  df1 <- df_wide
  df2 <- df_wide
  
  df1$is_A <- 1
  df1$is_B <- 0
  df1$AorB <- 1
  df1$pa_sub <- df1$pa_sub_A
  df1$pa_obj <- df1$pa_obj_A
  df1$aff <- df1$aff_A
  df1$reactance <- df1$reactance_A
  df1$anticipAff <- df1$anticipAff_A
  #df1$aff_reactance_pressure <- df1$aff_reactance_pressure_A
  #df1$aff_reactance_pushing <- df1$aff_reactance_pushing_A
  
  df2$is_A <- 0
  df2$is_B <- 1
  df2$AorB <- 2
  df2$pa_sub <- df2$pa_sub_B
  df2$pa_obj <- df2$pa_obj_B
  df2$aff <- df2$aff_B
  df2$reactance <- df2$reactance_B
  df2$anticipAff <- df2$anticipAff_B
  #df2$aff_reactance_pressure <- df2$aff_reactance_pressure_B
  #df2$aff_reactance_pushing <- df2$aff_reactance_pushing_B
  
  
  df_double <- rbind(df1, df2)

  
  # Double Stacking predictors
  # persuasion
  df_double$persuasion_self <- df_double$persuasion_A * df_double$is_A + df_double$persuasion_B * df_double$is_B
  df_double$persuasion_partner <- df_double$persuasion_A * df_double$is_B + df_double$persuasion_B * df_double$is_A
  
  df_double$persuasion_self_cw <- df_double$persuasion_cw_A * df_double$is_A + df_double$persuasion_cw_B * df_double$is_B
  df_double$persuasion_partner_cw <- df_double$persuasion_cw_A * df_double$is_B + df_double$persuasion_cw_B * df_double$is_A
  
  df_double$persuasion_self_cb <- df_double$persuasion_cb_A * df_double$is_A + df_double$persuasion_cb_B * df_double$is_B
  df_double$persuasion_partner_cb <- df_double$persuasion_cb_A * df_double$is_B + df_double$persuasion_cb_B * df_double$is_A
  
  #pressure
  df_double$pressure_self <- df_double$pressure_A * df_double$is_A + df_double$pressure_B * df_double$is_B
  df_double$pressure_partner <- df_double$pressure_A * df_double$is_B + df_double$pressure_B * df_double$is_A
  
  df_double$pressure_self_cw <- df_double$pressure_cw_A * df_double$is_A + df_double$pressure_cw_B * df_double$is_B
  df_double$pressure_partner_cw <- df_double$pressure_cw_A * df_double$is_B + df_double$pressure_cw_B * df_double$is_A
  
  df_double$pressure_self_cb <- df_double$pressure_cb_A * df_double$is_A + df_double$pressure_cb_B * df_double$is_B
  df_double$pressure_partner_cb <- df_double$pressure_cb_A * df_double$is_B + df_double$pressure_cb_B * df_double$is_A
  
  # pushing
  df_double$pushing_self <- df_double$pushing_A * df_double$is_A + df_double$pushing_B * df_double$is_B
  df_double$pushing_partner <- df_double$pushing_A * df_double$is_B + df_double$pushing_B * df_double$is_A
  
  df_double$pushing_self_cw <- df_double$pushing_cw_A * df_double$is_A + df_double$pushing_cw_B * df_double$is_B
  df_double$pushing_partner_cw <- df_double$pushing_cw_A * df_double$is_B + df_double$pushing_cw_B * df_double$is_A
  
  df_double$pushing_self_cb <- df_double$pushing_cb_A * df_double$is_A + df_double$pushing_cb_B * df_double$is_B
  df_double$pushing_partner_cb <- df_double$pushing_cb_A * df_double$is_B + df_double$pushing_cb_B * df_double$is_A
  
  # Weartime
  
  df_double$weartime_self <- df_double$weartime_A * df_double$is_A + df_double$weartime_B * df_double$is_B
  df_double$weartime_partner <- df_double$weartime_A * df_double$is_B + df_double$weartime_B * df_double$is_A
  
  df_double$weartime_self_cw <- df_double$weartime_cw_A * df_double$is_A + df_double$weartime_cw_B * df_double$is_B
  df_double$weartime_partner_cw <- df_double$weartime_cw_A * df_double$is_B + df_double$weartime_cw_B * df_double$is_A
  
  df_double$weartime_self_cb <- df_double$weartime_cb_A * df_double$is_A + df_double$weartime_cb_B * df_double$is_B
  df_double$weartime_partner_cb <- df_double$weartime_cb_A * df_double$is_B + df_double$weartime_cb_B * df_double$is_A
  
  # barriersiers
  df_double$barriers_self <- df_double$barriers_A * df_double$is_A + df_double$barriers_B * df_double$is_B
  df_double$barriers_partner <- df_double$barriers_A * df_double$is_B + df_double$barriers_B * df_double$is_A
  
  df_double$barriers_self_cw <- df_double$barriers_cw_A * df_double$is_A + df_double$barriers_cw_B * df_double$is_B
  df_double$barriers_partner_cw <- df_double$barriers_cw_A * df_double$is_B + df_double$barriers_cw_B * df_double$is_A
  
  df_double$barriers_self_cb <- df_double$barriers_cb_A * df_double$is_A + df_double$barriers_cb_B * df_double$is_B
  df_double$barriers_partner_cb <- df_double$barriers_cb_A * df_double$is_B + df_double$barriers_cb_B * df_double$is_A
  
  
  
  ## For Sensitivity Analyses necessary
  
  # Plans
  df_double$plan_self <- as.numeric(df_double$plan_A) * df_double$is_A + as.numeric(df_double$plan_B) * df_double$is_B
  df_double$plan_partner <- as.numeric(df_double$plan_A) * df_double$is_B + as.numeric(df_double$plan_B) * df_double$is_A
  
  # Support
  df_double$support_self <- df_double$support_A * df_double$is_A + df_double$support_B * df_double$is_B
  df_double$support_partner <- df_double$support_A * df_double$is_B + df_double$support_B * df_double$is_A
  
  df_double$support_self_cw <- df_double$support_cw_A * df_double$is_A + df_double$support_cw_B * df_double$is_B
  df_double$support_partner_cw <- df_double$support_cw_A * df_double$is_B + df_double$support_cw_B * df_double$is_A
  
  df_double$support_self_cb <- df_double$support_cb_A * df_double$is_A + df_double$support_cb_B * df_double$is_B
  df_double$support_partner_cb <- df_double$support_cb_A * df_double$is_B + df_double$support_cb_B * df_double$is_A
  
  # JITAIs
  df_double$got_JITAI_self <- df_double$got_JITAI_A * df_double$is_A + df_double$got_JITAI_B * df_double$is_B
  df_double$got_JITAI_partner <- df_double$got_JITAI_A * df_double$is_B + df_double$got_JITAI_B * df_double$is_A  
  
  
  df_double$pa_obj_log <- log(df_double$pa_obj)
      
  return(list(df_wide, df_double)) 
}


# And prepare original dataset
df_list <- prepare_df(df)

df_wide <- df_list[[1]]
df_double <- df_list[[2]]

Creating function to report measures

report_measures <- function(data, measures, ICC = TRUE) {
  
  if (ICC & any(sapply(data[,measures], is.factor))) {
    warning("Cannot add ICCs when Factors are in Dataframe.")
    ICC <- FALSE
  }
  
  measures_table <- as.data.frame(report::report_table(data[, measures]))

  measures_table <- measures_table %>%
    mutate(across(.cols = c('Mean', 'SD'),
                  .fns = ~format(round(.x, 2), nsmall = 2)))
  
  
  measures_table$Missing <- ifelse(
    is.na(measures_table$percentage_Missing),
    NA,
    paste0(round(measures_table$percentage_Missing), '%')
  )
  
  measures_table$Range <- ifelse(
    is.na(measures_table$Min) | is.na(measures_table$Max),
    NA,
    paste0(format(round(measures_table$Min, 2), nsmall = 2), '-', format(round(measures_table$Max, 2), nsmall = 2))
  )
    
  
  
  finished_df <- measures_table
  
  if (ICC) {
    cors <- wbCorr(data[,measures], data$coupleID)
    
    finished_df <- cbind(
      measures_table[, c('Variable', 'n_Obs', 'Missing', 'Mean', 'SD', 'Range')],
      ICC = format(round(get_icc(cors)$ICC, 2), nsmall = 2)
    )
  } else {
    measures_table$percentage_Obs <- ifelse(
      is.na(measures_table$percentage_Obs),
      NA,
      paste0(round(measures_table$percentage_Obs), '%')
    )
    
    finished_df <- measures_table[, c('Variable', 'Level', 'n_Obs', 'percentage_Obs', 'Missing', 'Mean', 'SD', 'Range')]
  }
  
  # Make them all Text for easier copying to excel. 
  finished_df[] <- lapply(finished_df, as.character)
  
  return(finished_df)
}

Sample Description

# mean of both relationship durations reported
df_wide$rel_duration_A <- df_wide$pre_rel_duration_y_A + df_wide$pre_rel_duration_m_A / 12 # in months
df_wide$rel_duration_B <- df_wide$pre_rel_duration_y_B + df_wide$pre_rel_duration_m_B / 12
df_wide$Relationship_Duration <- (df_wide$rel_duration_A + df_wide$rel_duration_B) / 2

# mean of both cohabiting duration reported
df_wide$hab_duration_A <- df_wide$pre_hab_duration_y_A + df_wide$pre_hab_duration_m_A / 12
df_wide$hab_duration_B <- df_wide$pre_hab_duration_y_B + df_wide$pre_hab_duration_m_B / 12 
df_wide$Cohabiting_Duration <- (df_wide$hab_duration_A + df_wide$hab_duration_B) / 2

#Gender of Each Partner
df_wide$Gender_A <- factor(df_wide$gender_A, levels = c(1,2,3), labels = c('Male','Female', 'Other'))
df_wide$Gender_B <- factor(df_wide$gender_B, levels = c(1,2,3), labels = c('Male','Female', 'Other'))
#Same sex couples or not
df_wide$Same_Sex <- factor(df_wide$Gender_A == df_wide$Gender_B, levels = c(FALSE, TRUE), labels = c('Same-Sex Couple', 'Mixed-Sex Couple'))

#Handedness
df_wide$Handedness_A <- factor(df_wide$pre_handedness_A, levels = c(0, 1, 2), c('Right','Left', 'Ambidextrous'))
df_wide$Handedness_B <- factor(df_wide$pre_handedness_B, levels = c(0, 1, 2), c('Right','Left', 'Ambidextrous'))



# Income (mean of both partner's report)

merge_income <- function(income1, income2) {
  merged_income <- numeric(length(income1))
  
  # Loop through each pair of incomes
  for (i in seq_along(income1)) {
    # Handle NA
    if (is.na(income1[i])) {
      merged_income[i] <- income2[i]
    } 
    
    else if (is.na(income2[i])) {
      merged_income[i] <- income1[i]
    }
    
    # if both are informative, take mean and round
    else if (income1[i] %in% 1:6 && income2[i] %in% 1:6) {
      merged_income[i] <- round((income1[i] + income2[i]) / 2)
    }
    
    # if one is informative and the other not, use the informative one. 
    else if (income1[i] %in% 1:6) {
      merged_income[i] <- income1[i]
    }
    else if (income2[i] %in% 1:6) {
      merged_income[i] <- income2[i]
    }
    
    # Now we only have cases, where both are either 70 or 99. We simply report "undisclosed". 
    else {
      merged_income[i] <- 99
    }
  }
  
  return(merged_income)
}


# Apply the function
numeric_income <- merge_income(df_wide$pre_income_1_A, df_wide$pre_income_1_B)

# Convert to factor
df_wide$Household_Income <- factor(numeric_income, levels = c(1,2,3,4,5,6,70,99), labels = c(
    "up to CHF 2'000.-", 
    "CHF 2'001.- to CHF 4'000.-",
    "CHF 4'001.- to CHF 6'000.-",
    "CHF 6'001.- to CHF 8'000.-",
    "CHF 8'001.- to CHF 10'000.-", 
    "above CHF 10'000.-", 
    "I don't know",
    "Undisclosed"
))





# Education for both separate
df_wide$Highest_Education_A <- factor(df_wide$pre_education_A, levels = c(1,2,3,4,5,6,7), labels = c(
  "(still) no school diploma",
  "compulsory education (9 years)",
  "vocational training (apprenticeship)",
  "Matura (university entrance qualification)",
  "Bachelor's degree", 
  "Master's degree",
  "Doctorate degree"
  )
)

df_wide$Highest_Education_B <- factor(df_wide$pre_education_B, levels = c(1,2,3,4,5,6,7), labels = c(
  "(still) no school diploma",
  "compulsory education (9 years)",
  "vocational training (apprenticeship)",
  "Matura (university entrance qualification)",
  "Bachelor's degree", 
  "Master's degree",
  "Doctorate degree"
  )
)



# Marital Status (married)
df_wide$Marital_Status_A <- df_wide$pre_mat_stat_A == 1
df_wide$Marital_status_B <- df_wide$pre_mat_stat_B == 1

# if at least one of them said they are married, we go for married
df_wide$Marital_Status <- factor((df_wide$Marital_Status_A + df_wide$Marital_status_B) > 0, levels = c(FALSE, TRUE), labels = c('Not Married', 'Married'))



# Children
df_wide$Children <- factor((df_wide$pre_child_option_A + df_wide$pre_child_option_B) > 0, levels = c(FALSE, TRUE), labels = c('Have Children', 'No Children'))
# number of chilren (ONLY FOR COUPLES THAT HAVE CHILDREN)
df_wide$nChildren_couples_with_children <- (df_wide$pre_child_nr_childs_A + df_wide$pre_child_nr_childs_B) / 2 
# number of children if every couple is considered. 
df_wide$nChildren_all_couples <- df_wide$nChildren_couples_with_children
df_wide$nChildren_all_couples[is.na(df_wide$nChildren_all_couples)] <- 0


# BMI (kg/m^2)
df_wide$pre_height_A <- df_wide$pre_height_A / 100 # to meters
df_wide$pre_height_B <- df_wide$pre_height_B / 100 # to meters

df_wide$BMI_A <- df_wide$pre_weight_A / (df_wide$pre_height_A^2)
df_wide$BMI_B <- df_wide$pre_weight_B / (df_wide$pre_height_B^2)



sample_measures <- c(
  "Gender_A", "Gender_B", "Same_Sex",
  "pre_age_A", "pre_age_B", "BMI_A", "BMI_B", "Handedness_A", "Handedness_B", # physical stats
  "Household_Income", "Highest_Education_A", "Highest_Education_B", # SES
  "Marital_Status", "Relationship_Duration", "Cohabiting_Duration",
  "Children", "nChildren_all_couples", "nChildren_couples_with_children" 
  )



sample_table <- report_measures(df_wide, sample_measures, ICC=F)


rownames(sample_table) <- NULL

openxlsx::write.xlsx(sample_table, 'Output/Sample.xlsx')

print_df(sample_table)
Variable Level n_Obs percentage_Obs Missing Mean SD Range
Gender_A Other 0 0% NA NA NA NA
Gender_A Female 825 39% NA NA NA NA
Gender_A Male 1265 61% NA NA NA NA
Gender_B Other 0 0% NA NA NA NA
Gender_B Male 825 39% NA NA NA NA
Gender_B Female 1265 61% NA NA NA NA
Same_Sex Mixed-Sex Couple 110 5% NA NA NA NA
Same_Sex Same-Sex Couple 1980 95% NA NA NA NA
pre_age_A NA 2090 NA 0% 34.05 11.19 19.00-58.00
pre_age_B NA 2090 NA 0% 33.97 10.72 19.00-60.00
BMI_A NA 2090 NA 0% 24.58 3.77 16.37-33.95
BMI_B NA 2090 NA 0% 25.29 4.33 17.32-33.95
Handedness_A Ambidextrous 55 3% NA NA NA NA
Handedness_A Left 330 16% NA NA NA NA
Handedness_A Right 1705 82% NA NA NA NA
Handedness_B Ambidextrous 0 0% NA NA NA NA
Handedness_B Left 275 13% NA NA NA NA
Handedness_B Right 1815 87% NA NA NA NA
Household_Income I don’t know 0 0% NA NA NA NA
Household_Income up to CHF 2’000.- 110 5% NA NA NA NA
Household_Income CHF 2’001.- to CHF 4’000.- 165 8% NA NA NA NA
Household_Income CHF 6’001.- to CHF 8’000.- 165 8% NA NA NA NA
Household_Income CHF 4’001.- to CHF 6’000.- 275 13% NA NA NA NA
Household_Income Undisclosed 275 13% NA NA NA NA
Household_Income CHF 8’001.- to CHF 10’000.- 440 21% NA NA NA NA
Household_Income above CHF 10’000.- 660 32% NA NA NA NA
Highest_Education_A (still) no school diploma 0 0% NA NA NA NA
Highest_Education_A compulsory education (9 years) 0 0% NA NA NA NA
Highest_Education_A Doctorate degree 0 0% NA NA NA NA
Highest_Education_A vocational training (apprenticeship) 330 16% NA NA NA NA
Highest_Education_A Master’s degree 495 24% NA NA NA NA
Highest_Education_A Matura (university entrance qualification) 605 29% NA NA NA NA
Highest_Education_A Bachelor’s degree 660 32% NA NA NA NA
Highest_Education_B (still) no school diploma 0 0% NA NA NA NA
Highest_Education_B compulsory education (9 years) 0 0% NA NA NA NA
Highest_Education_B Doctorate degree 0 0% NA NA NA NA
Highest_Education_B vocational training (apprenticeship) 110 5% NA NA NA NA
Highest_Education_B Bachelor’s degree 605 29% NA NA NA NA
Highest_Education_B Master’s degree 660 32% NA NA NA NA
Highest_Education_B Matura (university entrance qualification) 715 34% NA NA NA NA
Marital_Status Married 715 34% NA NA NA NA
Marital_Status Not Married 1375 66% NA NA NA NA
Relationship_Duration NA 2090 NA 0% 9.23 9.03 0.58-36.00
Cohabiting_Duration NA 2090 NA 0% 7.53 9.14 0.25-33.00
Children No Children 605 29% NA NA NA NA
Children Have Children 1485 71% NA NA NA NA
nChildren_all_couples NA 2090 NA 0% 0.47 0.85 0.00- 3.00
nChildren_couples_with_children NA 2090 NA 74% 1.80 0.60 1.00- 3.00

Descriptives of Measures

Main Measures

# Main Measures
main_constructs <- c("persuasion_A", "persuasion_B", 
                     "pressure_A", "pressure_B", 
                     "pushing_A", "pushing_B",
                     
                     "aff_A", "aff_B", 
                     "pa_sub_A", "pa_sub_B","pa_obj_A", "pa_obj_B", 
                     "anticipAff_A", "anticipAff_B", 
                     "reactance_A", "reactance_B"
                     )

main_descriptives <- report_measures(df_wide, main_constructs)

openxlsx::write.xlsx(main_descriptives, 'Output/DescriptivesMain.xlsx')

print_df(main_descriptives)
Variable n_Obs Missing Mean SD Range ICC
6 persuasion_A 2090 7% 0.43 1.09 0.00- 5.00 0.24
5 persuasion_B 2090 5% 0.41 1.08 0.00- 5.00 0.23
2 pressure_A 2090 7% 0.13 0.62 0.00- 5.00 0.63
1 pressure_B 2090 5% 0.12 0.62 0.00- 5.00 0.53
4 pushing_A 2090 7% 0.18 0.68 0.00- 5.00 0.09
3 pushing_B 2090 5% 0.15 0.64 0.00- 5.00 0.12
12 aff_A 2090 6% 4.89 1.13 1.00- 6.00 0.47
11 aff_B 2090 5% 4.77 1.14 1.00- 6.00 0.36
14 pa_sub_A 2090 6% 31.66 56.73 0.00-720.00 0.17
13 pa_sub_B 2090 5% 29.16 52.77 0.00-480.00 0.15
16 pa_obj_A 2090 12% 164.63 133.55 7.75-971.25 0.14
15 pa_obj_B 2090 11% 124.32 95.63 5.75-925.25 0.20
10 anticipAff_A 2090 6% 1.50 2.23 -5.00- 5.00 0.39
9 anticipAff_B 2090 5% 1.25 2.39 -5.00- 5.00 0.37
7 reactance_A 2090 80% 0.70 1.24 0.00- 4.00 0.52
8 reactance_B 2090 82% 0.90 1.39 0.00- 5.00 0.42

Including all Covariates

# With Covariates


# JITAIs as factor
df_wide$got_JITAI_A_factor <- factor(df_wide$got_JITAI_A, levels = c(0,1), labels = c("No JITAI", "Got a JITAI"))
df_wide$got_JITAI_B_factor <- factor(df_wide$got_JITAI_B, levels = c(0,1), labels = c("No JITAI", "Got a JITAI"))


# weekend as factor
df_wide$isWeekend_factor <- factor(df_wide$isWeekend_A, levels = c(0,1), labels = c("weekday", "weekend day"))


# Plan as a factor
df_wide$plan_A_factor <- factor(df_wide$plan_A, levels = c(0,1), labels = c("No Plan", "Yes Plan"))
df_wide$plan_B_factor <- factor(df_wide$plan_B, levels = c(0,1), labels = c("No Plan", "Yes Plan"))


# Study Group as a Factor
df_wide$studyGroup_factor <- factor(df_wide$studyGroup_A - 1, levels = c(0,1,2), labels = c("Allways Interventions", "First 3 weeks", "last 3 weeks"))

# Skilled support as factor
df_wide$skilled_support_factor <- factor(df_wide$skilled_support_A, levels = c(0,1), labels = c("Did not yet have SS", "Did already have SS"))



all_constructs <- c(
  main_constructs,
  "day",
  "weartime_A", "weartime_B",
  "barriers_A" , "barriers_B",
  "isWeekend_factor",
  "plan_A_factor", "plan_B_factor",
  "studyGroup_factor",
  "support_A", "support_B",
  "got_JITAI_A_factor", "got_JITAI_B_factor",
  "skilled_support_factor"
)


all_descriptives <- report_measures(df_wide, all_constructs, ICC = F)

openxlsx::write.xlsx(all_descriptives, 'Output/DescriptivesAll.xlsx')

print_df(all_descriptives)
Variable Level n_Obs percentage_Obs Missing Mean SD Range
29 persuasion_A NA 2090 NA 7% 0.43 1.09 0.00- 5.00
30 persuasion_B NA 2090 NA 5% 0.41 1.08 0.00- 5.00
31 pressure_A NA 2090 NA 7% 0.13 0.62 0.00- 5.00
32 pressure_B NA 2090 NA 5% 0.12 0.62 0.00- 5.00
33 pushing_A NA 2090 NA 7% 0.18 0.68 0.00- 5.00
34 pushing_B NA 2090 NA 5% 0.15 0.64 0.00- 5.00
18 aff_A NA 2090 NA 6% 4.89 1.13 1.00- 6.00
19 aff_B NA 2090 NA 5% 4.77 1.14 1.00- 6.00
27 pa_sub_A NA 2090 NA 6% 31.66 56.73 0.00- 720.00
28 pa_sub_B NA 2090 NA 5% 29.16 52.77 0.00- 480.00
25 pa_obj_A NA 2090 NA 12% 164.63 133.55 7.75- 971.25
26 pa_obj_B NA 2090 NA 11% 124.32 95.63 5.75- 925.25
20 anticipAff_A NA 2090 NA 6% 1.50 2.23 -5.00- 5.00
21 anticipAff_B NA 2090 NA 5% 1.25 2.39 -5.00- 5.00
35 reactance_A NA 2090 NA 80% 0.70 1.24 0.00- 4.00
36 reactance_B NA 2090 NA 82% 0.90 1.39 0.00- 5.00
24 day NA 2090 NA 0% 0.50 0.29 0.00- 1.00
39 weartime_A NA 2090 NA 2% 1036.85 374.95 0.00-1440.00
40 weartime_B NA 2090 NA 0% 1077.63 392.29 0.00-1440.00
22 barriers_A NA 2090 NA 54% 0.04 1.58 -3.29- 5.00
23 barriers_B NA 2090 NA 57% 0.10 1.70 -3.57- 5.00
6 isWeekend_factor weekend day 608 29% NA NA NA NA
14 isWeekend_factor weekday 1482 71% NA NA NA NA
2 plan_A_factor missing 136 7% NA NA NA NA
11 plan_A_factor Yes Plan 959 46% NA NA NA NA
12 plan_A_factor No Plan 995 48% NA NA NA NA
1 plan_B_factor missing 102 5% NA NA NA NA
10 plan_B_factor Yes Plan 901 43% NA NA NA NA
13 plan_B_factor No Plan 1087 52% NA NA NA NA
7 studyGroup_factor last 3 weeks 660 32% NA NA NA NA
8 studyGroup_factor Allways Interventions 715 34% NA NA NA NA
9 studyGroup_factor First 3 weeks 715 34% NA NA NA NA
37 support_A NA 2090 NA 7% 0.87 1.45 0.00- 5.00
38 support_B NA 2090 NA 5% 0.94 1.58 0.00- 5.00
3 got_JITAI_A_factor Got a JITAI 287 14% NA NA NA NA
17 got_JITAI_A_factor No JITAI 1803 86% NA NA NA NA
4 got_JITAI_B_factor Got a JITAI 298 14% NA NA NA NA
16 got_JITAI_B_factor No JITAI 1792 86% NA NA NA NA
5 skilled_support_factor Did not yet have SS 518 25% NA NA NA NA
15 skilled_support_factor Did already have SS 1572 75% NA NA NA NA

Correlations Main Constructs

cors <- wbCorr(df_wide[,c(main_constructs)], df_wide$coupleID, method = 'spearman')

main_cors <- summary(cors, 'wb')$merged_wb


openxlsx::write.xlsx(main_cors, 'Output/Correlations.xlsx')

print_df(main_cors, width = '7em')
persuasion_A persuasion_B pressure_A pressure_B pushing_A pushing_B aff_A aff_B pa_sub_A pa_sub_B pa_obj_A pa_obj_B anticipAff_A anticipAff_B reactance_A reactance_B
persuasion_A [0.24] 0.23*** 0.29*** 0.19*** 0.38*** 0.07** -0.02 0.03 0.14*** 0.11*** 0.01 0.00 0.02 0.03 -0.08 -0.01
persuasion_B 0.32 [0.23] 0.03 0.27*** 0.08*** 0.37*** -0.01 -0.01 0.20*** 0.20*** 0.07** 0.13*** 0.05* 0.04 0.01 0.02
pressure_A 0.46** 0.02 [0.63] 0.30*** 0.33*** -0.00 0.01 0.01 -0.01 -0.01 -0.03 -0.01 0.07** 0.06** 0.24*** 0.03
pressure_B 0.47** 0.34* 0.45** [0.53] 0.05* 0.23*** -0.02 -0.02 0.02 -0.04* -0.08** 0.06* 0.04 0.02 0.05 0.28***
pushing_A 0.55*** 0.23 0.48** 0.24 [0.09] 0.17*** -0.01 0.02 0.10*** 0.10*** 0.02 -0.02 0.02 0.02 0.09 -0.11
pushing_B 0.24 0.67*** 0.11 0.38* 0.35* [0.12] 0.03 0.02 0.17*** 0.13*** 0.07** 0.10*** 0.04 0.04 0.01 0.00
aff_A 0.35* 0.19 -0.03 -0.06 0.24 -0.03 [0.47] 0.19*** 0.19*** 0.11*** 0.13*** 0.12*** 0.30*** 0.10*** -0.02 0.00
aff_B 0.23 -0.04 0.01 -0.12 0.35* -0.01 0.41* [0.36] 0.10*** 0.20*** 0.10*** 0.18*** 0.08*** 0.36*** -0.02 -0.04
pa_sub_A 0.12 0.23 -0.27 -0.06 0.13 0.20 0.54*** 0.47** [0.17] 0.56*** 0.28*** 0.24*** 0.19*** 0.11*** -0.09 0.02
pa_sub_B 0.08 0.26 -0.16 -0.19 0.13 0.11 0.19 0.46** 0.70*** [0.15] 0.22*** 0.33*** 0.11*** 0.20*** 0.04 -0.05
pa_obj_A 0.03 0.22 0.02 -0.22 0.26 0.20 0.12 0.25 0.30 0.38* [0.14] 0.23*** 0.05* 0.03 0.01 0.03
pa_obj_B -0.06 0.18 0.04 -0.20 0.19 0.10 -0.00 0.40* 0.37* 0.54*** 0.76*** [0.20] 0.11*** 0.11*** 0.11* 0.13*
anticipAff_A 0.24 0.40* -0.35* -0.06 0.14 0.20 0.47** 0.38* 0.66*** 0.49** 0.40* 0.39* [0.39] 0.15*** -0.05 0.04
anticipAff_B 0.44** 0.12 0.07 0.25 0.39* 0.17 0.10 0.62*** 0.29 0.44** 0.11 0.28 0.41** [0.37] 0.00 -0.01
reactance_A 0.07 0.03 0.32* 0.07 0.14 0.03 -0.16 0.03 -0.13 -0.09 0.40* 0.37* 0.03 0.04 [0.52] 0.26***
reactance_B 0.17 0.30 -0.01 0.22 -0.16 0.29 -0.04 -0.15 -0.01 -0.03 0.07 0.06 0.11 0.02 0.01 [0.42]
plot(cors,'w')
## This may take a while...

#Modelling

Writing custion summary functions to facilitate reporting

summarize_brms <- function(model, short_version = FALSE, exponentiate = FALSE, model_rows_fixed = NULL, model_rows_random = NULL) {

  summ_og <- summary(model)
  summ <- summ_og$fixed
  rands <- summ_og$random[[1]][grep('sd\\(', rownames(summ_og$random[[1]])), ]
  
  ar1sterr <- summ_og$cor_pars
  sdresid <- summ_og$spec_pars
  
  rands <- rbind(rands, ar1sterr, sdresid)
  
  # Get rows that we need
  if (is.null(model_rows_fixed)) {
    model_rows_fixed <- rownames(summ)
  }
  if (is.null(model_rows_random)) {
    model_rows_random <- rownames(rands)
  }
  
  summ <- summ[model_rows_fixed, ]
  rands <- format(round(rands[model_rows_random, ], 2), nsmall = 2)
  
  Rhat <- format(round(summ$Rhat, 3), nsmall = 3)

  summ$Est.Error <- format(round(summ$Est.Error, 2), nsmall=2)
  summ$Bulk_ESS <- format(round(summ$Bulk_ESS, 2), nsmall=2)
  summ$Tail_ESS <- format(round(summ$Tail_ESS, 2), nsmall=2)
  
  #exponentiate if needed
  if (exponentiate) {
    summ$Estimate <- exp(summ$Estimate)
  }
    
  # Add stars based on Credible Interval
  significance <- ifelse(
    (summ$`l-95% CI` > 0 & summ$`u-95% CI` > 0) | 
    (summ$`l-95% CI` < 0 & summ$`u-95% CI` < 0), 
    '*', 
    '')
  
  # exponentiate CI
  if (exponentiate) {
    summ$`u-95% CI`<- exp(summ$`u-95% CI`)
    summ$`l-95% CI` <- exp(summ$`l-95% CI`)
  }
      

  
  # Round CI
  summ$`l-95% CI` <- format(round(summ$`l-95% CI`, 2), nsmall = 2)
  summ$`u-95% CI` <- format(round(summ$`u-95% CI`, 2), nsmall = 2)
  
  
  summ$Estimate <- ifelse( 
    is.na(summ$Estimate), 
    NA, 
    paste0(
      format(round(summ$Estimate,2), nsmall=2), 
      significance
    )
  )
  
  
  # Rename "Estimate" column correctly
  if (exponentiate) {
    if ('bernoulli' %in% model$family) {
      correct_name <- 'OR'
    } else if ('negbinomial' %in% model$family) {
      correct_name <- 'IRR'
    } else {
      correct_name <- 'exp(Est.)'
      warning(
        "Coefficients were exponentiated. Double check if this was intended."
        )
    }
    
  } else {
    correct_name <- 'b'
  }
  
  colnames(summ)[1] <- correct_name
  colnames(rands)[1] <- correct_name
  
  
  # Add rhat rounded to 3 digits back
  summ$Rhat <- Rhat
  
  # Add back all rownames
  rownames(summ) <- model_rows_fixed
  rownames(rands) <- model_rows_random
  
  # Add random effects to fixed effects df. 
  fullrep <- rbind(summ, rands)
  

  
  # We are finished, if we want the full table
  if (!short_version) {
      return(fullrep[ , c(1, 3, 4, 5, 6, 7)])
  }
  
  
  fullrep$`95% CI` <-  ifelse(
    is.na(fullrep[correct_name]) | fullrep$`u-95% CI` == "NA",
    NA,
    paste0("[", fullrep$`l-95% CI`, ", ", fullrep$`u-95% CI`, "]")
  )

  return(fullrep[ , c(1, 8)])
}


report_side_by_side <- function(..., model_rows_random, model_rows_fixed) {
  models <- list(...)
  model_names <- sapply(substitute(list(...))[-1], deparse)
  names(models) <- model_names

  side_by_side <- NULL
  for (i in seq_along(models)) {
    model <- models[[i]]
    model_name <- model_names[i]
    print(model_name)
    if ('bernoulli' %in% model$family | 'negbinomial' %in% model$family | grepl('log', model_name)) {
      exponentiate <- TRUE
    } else {
      exponentiate <- FALSE
    }
    model_summary <- summarize_brms(
      model, short_version = TRUE, 
      exponentiate = exponentiate, 
      model_rows_random = model_rows_random,
      model_rows_fixed = model_rows_fixed
      )
    
    colnames(model_summary) <- paste(colnames(model_summary), model_name)
    
    if (is.null(side_by_side)) {
      side_by_side <- model_summary
    } else {
      side_by_side <- cbind(side_by_side, model_summary)
    }
  }
  return(side_by_side)
}

MODELLING Persuasion

# For indistinguishable Dyads
model_rows_fixed_persuasion_ind <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'day', 
    'weartime_self_cw', 
    'barriers_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'weartime_self_cb',
    'barriers_self_cb'
  )
  
model_rows_random_persuasion_ind <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)


# For APIMs

model_rows_fixed_persuasion_apim <- c(
    'is_A',
    'is_B',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'is_A:persuasion_self_cw',
    'is_B:persuasion_self_cw', 
    'is_A:persuasion_partner_cw', 
    'is_B:persuasion_partner_cw',
    'is_A:day', 
    'is_B:day', 
    'is_A:weartime_self_cw', 
    'is_B:weartime_self_cw', 
    'is_A:barriers_self_cw',
    'is_B:barriers_self_cw',
    
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'is_A:persuasion_self_cb',
    'is_B:persuasion_self_cb', 
    'is_A:persuasion_partner_cb', 
    'is_B:persuasion_partner_cb',
    'is_A:weartime_self_cb', 
    'is_B:weartime_self_cb', 
    'is_A:barriers_self_cb',
    'is_B:barriers_self_cb'
  )
  

model_rows_random_persuasion_apim <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(is_A)',
  'sd(is_B)',
  'sd(is_A:persuasion_self_cw)',
  'sd(is_B:persuasion_self_cw)', 
  'sd(is_A:persuasion_partner_cw)', 
  'sd(is_B:persuasion_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

Subjective MVPA ~ persuasion

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 100)

Modelling using the gaussian family fails. Due to the many zeros, transformations won’t help estimating the models. We employ the negative binomial family.

Indistinguishable

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day + 
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_sub_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/persuasion_pa_sub_ind"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(persuasion_pa_sub_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(persuasion_pa_sub_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_pa_sub_ind)
## Warning: Found 1 observations with a pareto_k > 0.7 in
## model 'persuasion_pa_sub_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
print(loo_ind)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12076.4 177.6
## p_loo        23.7   2.4
## looic     24152.8 355.1
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 2.6]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3741  100.0%  438     
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.
plot(persuasion_pa_sub_ind, ask = FALSE)

summarize_brms(
  persuasion_pa_sub_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 27.55* 20.81 36.69 1.002 2331.38 2837.67
Within-Person Effects
persuasion_self_cw 1.24* 1.11 1.40 1.002 5869.29 3447.36
persuasion_partner_cw 1.25* 1.13 1.41 1.000 6134.98 3199.10
day 0.78 0.57 1.09 1.000 8469.49 3183.32
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb 1.06 0.76 1.49 1.002 5999.27 3231.88
persuasion_partner_cb 0.79 0.57 1.10 0.999 6336.80 3474.68
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.67 0.50 0.89 1.00 1445.36 2613.93
sd(persuasion_self_cw) 0.12 0.01 0.26 1.00 2432.34 2261.02
sd(persuasion_partner_cw) 0.10 0.00 0.24 1.00 2142.57 1879.26
Additional Parameters
ar[1] 0.02 -0.92 0.93 1.00 6362.35 2484.17
nu NA NA NA NA NA NA
shape 0.14 0.13 0.14 1.01 8960.37 2726.57
sderr 0.05 0.00 0.13 1.00 3223.60 2619.73
sigma NA NA NA NA NA NA

APIM

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_sub_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/persuasion_pa_sub_apim"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(persuasion_pa_sub_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(persuasion_pa_sub_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_pa_sub_apim)
## Warning: Found 2 observations with a pareto_k > 0.7 in
## model 'persuasion_pa_sub_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12074.7 177.6
## p_loo        37.4   3.3
## looic     24149.3 355.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3740  99.9%   341     
##    (0.7, 1]   (bad)         2   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(persuasion_pa_sub_apim, ask = FALSE)

summarize_brms(
  persuasion_pa_sub_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 30.51* 21.31 43.56 1.000 1890.46 2379.54
is_B 25.57* 18.34 35.70 1.001 2368.85 2596.26
Within-Person Effects
is_A:persuasion_self_cw 1.25* 1.07 1.48 1.000 5215.74 2929.61
is_B:persuasion_self_cw 1.25* 1.08 1.46 1.000 5386.14 2881.34
is_A:persuasion_partner_cw 1.24* 1.07 1.45 1.002 6066.25 2701.75
is_B:persuasion_partner_cw 1.28* 1.10 1.55 1.002 4003.64 2802.29
is_A:day 0.66 0.43 1.04 1.000 5545.58 2746.90
is_B:day 0.87 0.55 1.37 1.000 4771.13 2959.81
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb 0.88 0.43 1.85 1.002 1448.62 2263.63
is_B:persuasion_self_cb 1.45 0.72 2.95 1.001 1880.36 2501.68
is_A:persuasion_partner_cb 1.14 0.50 2.66 1.002 1490.78 2376.29
is_B:persuasion_partner_cb 0.60 0.32 1.14 1.001 1952.12 2305.28
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.78 0.56 1.05 1.00 1793.60 2162.97
sd(is_B) 0.63 0.44 0.86 1.00 2247.82 2610.67
sd(is_A:persuasion_self_cw) 0.15 0.01 0.33 1.00 2014.51 2042.92
sd(is_B:persuasion_self_cw) 0.11 0.01 0.30 1.00 1795.81 2219.22
sd(is_A:persuasion_partner_cw) 0.08 0.00 0.24 1.00 2678.70 2128.12
sd(is_B:persuasion_partner_cw) 0.18 0.02 0.39 1.00 2136.82 2165.95
Additional Parameters
ar[1] 0.03 -0.93 0.93 1.00 5219.11 2841.53
nu NA NA NA NA NA NA
shape 0.14 0.13 0.15 1.00 7643.78 2700.45
sderr 0.05 0.00 0.13 1.00 2500.90 1998.53
sigma NA NA NA NA NA NA

Compare models

loo_compare(loo_ind, loo_apim)
##                        elpd_diff se_diff
## persuasion_pa_sub_apim  0.0       0.0   
## persuasion_pa_sub_ind  -1.7       3.3
bayes_R2(persuasion_pa_sub_apim)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2541458 0.06157753 0.1575417 0.4013804
bayes_R2(persuasion_pa_sub_ind)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.1832215 0.04449012 0.1123851 0.2813885

Device Based MVPA ~ persuasion

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

We tried negative binomial here as well for consistency, but the model did not converge. Poisson also did not work. As we have no zeros in this distribution, we log transform.

Indistinguishable

formula <- bf(
  pa_obj_log ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_obj_log_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/persuasion_pa_obj_log_ind"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(persuasion_pa_obj_log_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(persuasion_pa_obj_log_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_pa_obj_log_ind)
print(loo_ind)
## 
## Computed from 4000 by 3305 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2796.8  55.8
## p_loo        68.2   3.0
## looic      5593.6 111.7
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.6]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_pa_obj_log_ind, ask = FALSE)

summarize_brms(
  persuasion_pa_obj_log_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
## Warning in summarize_brms(persuasion_pa_obj_log_ind,
## model_rows_fixed = model_rows_fixed_persuasion_ind, :
## Coefficients were exponentiated. Double check if this
## was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 117.72* 105.27 131.58 1.008 574.58 1215.75
Within-Person Effects
persuasion_self_cw 1.03 1.00 1.06 1.001 2062.72 2504.79
persuasion_partner_cw 1.02 0.99 1.05 1.001 2767.45 2970.09
day 0.97 0.89 1.06 1.001 6813.07 2674.23
weartime_self_cw 1.00* 1.00 1.00 1.003 4085.95 2541.69
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb 1.06 0.93 1.20 1.002 2128.14 2523.86
persuasion_partner_cb 1.01 0.89 1.14 1.000 2217.26 2933.59
weartime_self_cb 1.00 1.00 1.00 1.001 2999.76 3219.09
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.29 0.23 0.38 1.00 826.63 1626.99
sd(persuasion_self_cw) 0.05 0.03 0.08 1.00 1600.29 1324.44
sd(persuasion_partner_cw) 0.05 0.02 0.08 1.00 1874.29 1513.73
Additional Parameters
ar[1] 0.30 0.26 0.33 1.00 5355.39 3187.70
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.56 0.54 0.57 1.00 6994.95 3088.34

apim

formula <- bf(
  pa_obj_log ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_obj_log_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/persuasion_pa_obj_log_apim"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(persuasion_pa_obj_log_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(persuasion_pa_obj_log_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_pa_obj_log_apim)
print(loo_apim)
## 
## Computed from 4000 by 3305 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2722.4  56.7
## p_loo       120.8   5.6
## looic      5444.8 113.4
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.3]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_pa_obj_log_apim, ask = FALSE)

summarize_brms(
  persuasion_pa_obj_log_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
## Warning in summarize_brms(persuasion_pa_obj_log_apim,
## model_rows_fixed = model_rows_fixed_persuasion_apim, :
## Coefficients were exponentiated. Double check if this
## was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 137.00* 121.54 153.74 1.000 1146.59 1833.07
is_B 101.65* 89.08 116.52 1.002 1387.52 2279.23
Within-Person Effects
is_A:persuasion_self_cw 1.00 0.97 1.04 1.001 3317.88 2987.03
is_B:persuasion_self_cw 1.05* 1.01 1.09 1.000 3417.72 3135.24
is_A:persuasion_partner_cw 1.03 0.99 1.07 1.000 3443.29 2766.30
is_B:persuasion_partner_cw 1.02 0.98 1.05 1.002 3558.63 2919.08
is_A:day 0.94 0.83 1.06 1.007 6249.65 3283.95
is_B:day 1.00 0.89 1.12 1.002 6139.29 3073.27
is_A:weartime_self_cw 1.00* 1.00 1.00 1.000 4398.93 2468.88
is_B:weartime_self_cw 1.00* 1.00 1.00 1.002 4207.09 2774.05
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb 0.88 0.68 1.13 1.001 1051.16 2059.39
is_B:persuasion_self_cb 1.28 0.90 1.79 1.001 1311.70 1876.12
is_A:persuasion_partner_cb 1.29 0.97 1.72 1.003 1019.86 1565.15
is_B:persuasion_partner_cb 0.88 0.66 1.18 1.001 1249.57 2068.82
is_A:weartime_self_cb 1.00 1.00 1.00 1.002 1523.84 2272.45
is_B:weartime_self_cb 1.00 1.00 1.00 1.000 2170.84 3073.18
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.28 0.21 0.37 1.00 1474.91 2148.46
sd(is_B) 0.34 0.26 0.44 1.01 1058.25 2392.37
sd(is_A:persuasion_self_cw) 0.06 0.01 0.10 1.00 1395.30 1242.92
sd(is_B:persuasion_self_cw) 0.06 0.02 0.11 1.00 1660.07 1144.72
sd(is_A:persuasion_partner_cw) 0.06 0.01 0.12 1.00 1071.24 1157.70
sd(is_B:persuasion_partner_cw) 0.05 0.01 0.10 1.00 1334.76 1193.59
Additional Parameters
ar[1] 0.23 0.19 0.26 1.00 5693.73 2771.04
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.54 0.53 0.55 1.00 6557.80 2705.73

Compare Models

loo_compare(loo_ind, loo_apim)
##                            elpd_diff se_diff
## persuasion_pa_obj_log_apim   0.0       0.0  
## persuasion_pa_obj_log_ind  -74.4      14.1
bayes_R2(persuasion_pa_obj_log_apim)
##     Estimate   Est.Error      Q2.5    Q97.5
## R2 0.1872719 0.006202424 0.1748286 0.199127
bayes_R2(persuasion_pa_obj_log_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.1681875 0.006551877 0.1555223 0.1810579

Affect ~ persuasion

range(df_double$aff, na.rm = T) 
## [1] 1 6
hist(df_double$aff, breaks = 15)

Indistinguishable

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_aff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/persuasion_aff_ind"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(persuasion_aff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(persuasion_aff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_aff_ind)
print(loo_ind)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4823.6  63.5
## p_loo        62.9   2.8
## looic      9647.2 126.9
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_aff_ind, ask = FALSE)

summarize_brms(
  persuasion_aff_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.72* 4.52 4.93 1.002 625.85 1554.29
Within-Person Effects
persuasion_self_cw 0.00 -0.03 0.03 1.000 4191.05 2934.22
persuasion_partner_cw 0.03 -0.01 0.07 1.001 3794.74 3145.50
day 0.21* 0.05 0.38 1.001 7024.33 2584.35
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb 0.13 -0.09 0.34 1.001 2252.92 2739.00
persuasion_partner_cb 0.08 -0.14 0.30 1.000 2369.12 2921.80
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.60 0.47 0.78 1.00 945.11 1813.70
sd(persuasion_self_cw) 0.03 0.00 0.08 1.00 1258.72 1240.37
sd(persuasion_partner_cw) 0.07 0.01 0.12 1.00 1165.99 1341.53
Additional Parameters
ar[1] 0.45 0.42 0.48 1.00 5361.23 3387.10
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.87 0.85 0.89 1.00 6099.35 3060.03

APIM

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_aff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/persuasion_aff_apim"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(persuasion_aff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(persuasion_aff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_aff_apim)
## Warning: Found 1 observations with a pareto_k > 0.7 in
## model 'persuasion_aff_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4684.2  65.5
## p_loo       119.8   5.4
## looic      9368.5 131.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3741  100.0%  527     
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.
plot(persuasion_aff_apim, ask = FALSE)

summarize_brms(
  persuasion_aff_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 4.84* 4.55 5.14 1.016 319.15 625.53
is_B 4.56* 4.32 4.81 1.003 631.82 1087.51
Within-Person Effects
is_A:persuasion_self_cw -0.01 -0.05 0.04 1.000 2961.27 3077.12
is_B:persuasion_self_cw 0.01 -0.05 0.06 1.001 2577.34 2764.01
is_A:persuasion_partner_cw 0.01 -0.04 0.06 1.003 3559.99 2911.73
is_B:persuasion_partner_cw 0.05 -0.01 0.10 1.001 2008.47 2632.16
is_A:day 0.08 -0.10 0.26 1.001 4551.98 2857.39
is_B:day 0.40* 0.20 0.58 1.003 4114.48 3139.34
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb 0.54 -0.15 1.22 1.004 469.33 996.31
is_B:persuasion_self_cb -0.50 -1.22 0.20 1.010 517.28 993.26
is_A:persuasion_partner_cb -0.40 -1.24 0.45 1.007 346.29 694.75
is_B:persuasion_partner_cb 0.51 -0.08 1.11 1.008 494.53 1005.98
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.79 0.60 1.03 1.01 797.90 1121.53
sd(is_B) 0.70 0.54 0.91 1.00 959.13 1692.08
sd(is_A:persuasion_self_cw) 0.05 0.00 0.13 1.00 833.26 1366.44
sd(is_B:persuasion_self_cw) 0.08 0.02 0.16 1.01 985.60 680.38
sd(is_A:persuasion_partner_cw) 0.06 0.00 0.13 1.00 1032.15 1659.39
sd(is_B:persuasion_partner_cw) 0.10 0.03 0.17 1.00 1038.19 785.70
Additional Parameters
ar[1] 0.33 0.29 0.36 1.00 3881.52 3131.03
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.83 0.81 0.85 1.00 4082.76 2917.19

Compare Models

loo_compare(loo_ind, loo_apim)
##                     elpd_diff se_diff
## persuasion_aff_apim    0.0       0.0 
## persuasion_aff_ind  -139.4      19.9
bayes_R2(persuasion_aff_apim)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2358568 0.004479439 0.2268888 0.2445024
bayes_R2(persuasion_aff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2240658 0.004720843 0.2147554 0.2333211

Anticipatory Affect ~ persuasion

range(df_double$anticipAff, na.rm = T) 
## [1] -5  5
hist(df_double$anticipAff, breaks = 10)

Indistinguishable

formula <- bf(
  anticipAff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=-5, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_anticipAff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/persuasion_anticipAff_ind"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(persuasion_anticipAff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(persuasion_anticipAff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_anticipAff_ind)
print(loo_ind)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7522.4  56.4
## p_loo        62.4   2.6
## looic     15044.9 112.9
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.7]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_anticipAff_ind, ask = FALSE)

summarize_brms(
  persuasion_anticipAff_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.48* 1.03 1.91 1.003 680.82 1183.55
Within-Person Effects
persuasion_self_cw 0.06 -0.03 0.14 1.001 2359.76 2281.74
persuasion_partner_cw 0.08* 0.02 0.14 1.000 4177.68 2400.13
day -0.18 -0.52 0.14 1.001 5043.71 3161.45
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb -0.31 -0.74 0.15 1.001 2032.85 2383.06
persuasion_partner_cb 1.01* 0.57 1.44 1.000 2003.70 2647.24
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.22 0.95 1.58 1.01 779.45 1610.71
sd(persuasion_self_cw) 0.17 0.08 0.27 1.00 1194.87 1670.16
sd(persuasion_partner_cw) 0.03 0.00 0.10 1.00 2181.81 1809.49
Additional Parameters
ar[1] 0.41 0.38 0.44 1.00 4576.52 2904.04
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.79 1.75 1.83 1.00 5012.64 2726.37

APIM

formula <- bf(
  anticipAff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_anticipAff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 5000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/persuasion_anticipAff_apim"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(persuasion_anticipAff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(persuasion_anticipAff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_anticipAff_apim)
print(loo_apim)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7435.8  56.6
## p_loo       109.3   4.2
## looic     14871.6 113.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.0]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_anticipAff_apim, ask = FALSE)

summarize_brms(
  persuasion_anticipAff_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 1.53* 1.03 2.02 1.014 587.58 956.01
is_B 1.42* 0.91 1.92 1.002 684.20 1395.61
Within-Person Effects
is_A:persuasion_self_cw 0.07 -0.03 0.18 1.001 3166.66 2642.78
is_B:persuasion_self_cw 0.04 -0.07 0.16 1.001 2224.67 2949.46
is_A:persuasion_partner_cw 0.07 -0.02 0.16 1.001 4626.46 2986.92
is_B:persuasion_partner_cw 0.08 -0.01 0.16 1.002 4316.70 3153.17
is_A:day -0.01 -0.39 0.39 1.001 5160.00 2801.20
is_B:day -0.35 -0.74 0.04 1.001 4584.55 3332.61
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb -0.38 -1.53 0.83 1.003 691.02 1398.90
is_B:persuasion_self_cb -0.41 -1.96 0.98 1.005 701.57 1510.27
is_A:persuasion_partner_cb 1.16 -0.16 2.55 1.001 633.55 1116.65
is_B:persuasion_partner_cb 0.92 -0.30 2.27 1.006 710.45 1744.78
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 1.40 1.09 1.83 1.00 670.15 1221.52
sd(is_B) 1.46 1.14 1.87 1.01 912.46 1854.03
sd(is_A:persuasion_self_cw) 0.14 0.01 0.28 1.00 1035.98 1259.52
sd(is_B:persuasion_self_cw) 0.21 0.08 0.35 1.00 1460.96 1285.15
sd(is_A:persuasion_partner_cw) 0.07 0.00 0.19 1.00 1636.37 1541.32
sd(is_B:persuasion_partner_cw) 0.05 0.00 0.16 1.00 1858.63 1743.17
Additional Parameters
ar[1] 0.33 0.30 0.37 1.00 4490.48 2819.45
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.74 1.70 1.78 1.00 6484.89 2619.98

Compare Models

loo_compare(loo_ind, loo_apim)
##                            elpd_diff se_diff
## persuasion_anticipAff_apim   0.0       0.0  
## persuasion_anticipAff_ind  -86.7      14.7
bayes_R2(persuasion_anticipAff_apim)
##     Estimate   Est.Error     Q2.5     Q97.5
## R2 0.2186067 0.004412307 0.209873 0.2271955
bayes_R2(persuasion_anticipAff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2112572 0.004741973 0.2020313 0.2206235

reactance ~ persuasion

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 6)

Indistinguishable

formula <- bf(
  reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_reactance_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/persuasion_reactance_ind"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(persuasion_reactance_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(persuasion_reactance_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_reactance_ind)
## Warning: Found 1 observations with a pareto_k > 0.7 in
## model 'persuasion_reactance_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
print(loo_ind)
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1097.6 34.0
## p_loo        51.2  5.3
## looic      2195.2 68.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.8]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     755   99.9%   215     
##    (0.7, 1]   (bad)        1    0.1%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(persuasion_reactance_ind, ask = FALSE)

summarize_brms(
  persuasion_reactance_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.53* 0.31 0.75 1.001 3837.80 2918.10
Within-Person Effects
persuasion_self_cw -0.03 -0.10 0.04 1.001 4054.41 3432.26
persuasion_partner_cw 0.01 -0.05 0.08 1.001 5647.22 3422.90
day 0.17 -0.12 0.46 1.001 7488.77 2985.47
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb 0.48* 0.26 0.71 1.000 5029.33 2530.32
persuasion_partner_cb 0.34* 0.12 0.56 1.001 6454.69 2974.56
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.41 0.26 0.60 1.00 1792.58 2582.37
sd(persuasion_self_cw) 0.12 0.02 0.23 1.01 770.13 1084.24
sd(persuasion_partner_cw) 0.06 0.00 0.14 1.00 1974.15 2584.17
Additional Parameters
ar[1] 0.06 -0.02 0.14 1.00 5726.54 3042.33
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.00 0.95 1.05 1.00 4152.46 2245.30

APIM

formula <- bf(
  reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 

  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_reactance_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/persuasion_reactance_apim"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(persuasion_reactance_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(persuasion_reactance_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_reactance_apim)
## Warning: Found 3 observations with a pareto_k > 0.7 in
## model 'persuasion_reactance_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1080.9 33.0
## p_loo        83.9  7.9
## looic      2161.8 66.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.6]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     753   99.6%   169     
##    (0.7, 1]   (bad)        3    0.4%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(persuasion_reactance_apim, ask = FALSE)

summarize_brms(
  persuasion_reactance_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 0.63* 0.34 0.92 1.004 1553.91 2419.08
is_B 0.44* 0.18 0.72 1.001 2545.34 2746.90
Within-Person Effects
is_A:persuasion_self_cw -0.07 -0.19 0.04 1.004 1579.07 2194.33
is_B:persuasion_self_cw 0.00 -0.08 0.09 1.000 2672.03 2970.44
is_A:persuasion_partner_cw 0.09 0.00 0.18 1.002 3245.89 3044.42
is_B:persuasion_partner_cw -0.06 -0.15 0.04 1.001 3599.14 2866.75
is_A:day -0.04 -0.42 0.31 1.001 3371.65 2271.19
is_B:day 0.32 -0.08 0.71 1.002 3972.40 3106.95
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb 0.51 -0.06 1.10 1.003 903.67 1449.95
is_B:persuasion_self_cb 0.56* 0.05 1.05 1.002 1569.23 2313.50
is_A:persuasion_partner_cb 0.37 -0.26 1.01 1.002 1128.94 1792.12
is_B:persuasion_partner_cb 0.29 -0.14 0.73 1.002 1191.99 1575.85
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.59 0.38 0.86 1.00 1152.44 1880.48
sd(is_B) 0.41 0.21 0.66 1.00 1233.58 2055.63
sd(is_A:persuasion_self_cw) 0.23 0.11 0.37 1.00 840.20 976.89
sd(is_B:persuasion_self_cw) 0.09 0.01 0.21 1.01 795.45 1859.45
sd(is_A:persuasion_partner_cw) 0.09 0.01 0.20 1.00 1580.36 1689.74
sd(is_B:persuasion_partner_cw) 0.09 0.01 0.21 1.00 1553.30 1922.96
Additional Parameters
ar[1] -0.03 -0.12 0.06 1.00 2565.53 2596.32
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.95 0.90 1.01 1.00 2581.88 2938.75

Compare Models

loo_compare(loo_ind, loo_apim)
##                           elpd_diff se_diff
## persuasion_reactance_apim   0.0       0.0  
## persuasion_reactance_ind  -16.7       7.1
bayes_R2(persuasion_reactance_apim)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2358401 0.01024904 0.2147279 0.2547377
bayes_R2(persuasion_reactance_ind)
##     Estimate  Est.Error     Q2.5     Q97.5
## R2 0.2185654 0.01072507 0.197132 0.2390404

Report models for Persuasion side by side

Indistinguishable Models

persuasion_all_ind_models <- report_side_by_side(
  persuasion_pa_sub_ind,
  persuasion_pa_obj_log_ind,
  persuasion_aff_ind,
  persuasion_anticipAff_ind,
  persuasion_reactance_ind,
  
  model_rows_random = model_rows_random_persuasion_ind,
  model_rows_fixed = model_rows_fixed_persuasion_ind
) 
## [1] "persuasion_pa_sub_ind"
## [1] "persuasion_pa_obj_log_ind"
## Warning in summarize_brms(model, short_version = TRUE,
## exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "persuasion_aff_ind"
## [1] "persuasion_anticipAff_ind"
## [1] "persuasion_reactance_ind"
# pretty printing
persuasion_all_ind_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_ind()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR persuasion_pa_sub_ind 95% CI persuasion_pa_sub_ind exp(Est.) persuasion_pa_obj_log_ind 95% CI persuasion_pa_obj_log_ind b persuasion_aff_ind 95% CI persuasion_aff_ind b persuasion_anticipAff_ind 95% CI persuasion_anticipAff_ind b persuasion_reactance_ind 95% CI persuasion_reactance_ind
Intercept 27.55* [20.81, 36.69] 117.72* [105.27, 131.58] 4.72* [ 4.52, 4.93] 1.48* [ 1.03, 1.91] 0.53* [ 0.31, 0.75]
Within-Person Effects
persuasion_self_cw 1.24* [ 1.11, 1.40] 1.03 [ 1.00, 1.06] 0.00 [-0.03, 0.03] 0.06 [-0.03, 0.14] -0.03 [-0.10, 0.04]
persuasion_partner_cw 1.25* [ 1.13, 1.41] 1.02 [ 0.99, 1.05] 0.03 [-0.01, 0.07] 0.08* [ 0.02, 0.14] 0.01 [-0.05, 0.08]
day 0.78 [ 0.57, 1.09] 0.97 [ 0.89, 1.06] 0.21* [ 0.05, 0.38] -0.18 [-0.52, 0.14] 0.17 [-0.12, 0.46]
weartime_self_cw NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb 1.06 [ 0.76, 1.49] 1.06 [ 0.93, 1.20] 0.13 [-0.09, 0.34] -0.31 [-0.74, 0.15] 0.48* [ 0.26, 0.71]
persuasion_partner_cb 0.79 [ 0.57, 1.10] 1.01 [ 0.89, 1.14] 0.08 [-0.14, 0.30] 1.01* [ 0.57, 1.44] 0.34* [ 0.12, 0.56]
weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.67 [ 0.50, 0.89] 0.29 [0.23, 0.38] 0.60 [0.47, 0.78] 1.22 [0.95, 1.58] 0.41 [ 0.26, 0.60]
sd(persuasion_self_cw) 0.12 [ 0.01, 0.26] 0.05 [0.03, 0.08] 0.03 [0.00, 0.08] 0.17 [0.08, 0.27] 0.12 [ 0.02, 0.23]
sd(persuasion_partner_cw) 0.10 [ 0.00, 0.24] 0.05 [0.02, 0.08] 0.07 [0.01, 0.12] 0.03 [0.00, 0.10] 0.06 [ 0.00, 0.14]
Additional Parameters
ar[1] 0.02 [-0.92, 0.93] 0.30 [0.26, 0.33] 0.45 [0.42, 0.48] 0.41 [0.38, 0.44] 0.06 [-0.02, 0.14]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.14 [ 0.13, 0.14] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.05 [ 0.00, 0.13] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.56 [0.54, 0.57] 0.87 [0.85, 0.89] 1.79 [1.75, 1.83] 1.00 [ 0.95, 1.05]
# prepare export to excel
all_models <- list()
all_models$Persuasion_ind <- prepare_for_export(persuasion_all_ind_models)

APIMs

persuasion_all_apim_models <- report_side_by_side(
  persuasion_pa_sub_apim,
  persuasion_pa_obj_log_apim,
  persuasion_aff_apim,
  persuasion_anticipAff_apim,
  persuasion_reactance_apim,
  
  model_rows_random = model_rows_random_persuasion_apim,
  model_rows_fixed = model_rows_fixed_persuasion_apim
) 
## [1] "persuasion_pa_sub_apim"
## [1] "persuasion_pa_obj_log_apim"
## Warning in summarize_brms(model, short_version = TRUE,
## exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "persuasion_aff_apim"
## [1] "persuasion_anticipAff_apim"
## [1] "persuasion_reactance_apim"
# pretty printing
persuasion_all_apim_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_apim()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR persuasion_pa_sub_apim 95% CI persuasion_pa_sub_apim exp(Est.) persuasion_pa_obj_log_apim 95% CI persuasion_pa_obj_log_apim b persuasion_aff_apim 95% CI persuasion_aff_apim b persuasion_anticipAff_apim 95% CI persuasion_anticipAff_apim b persuasion_reactance_apim 95% CI persuasion_reactance_apim
is_A 30.51* [21.31, 43.56] 137.00* [121.54, 153.74] 4.84* [ 4.55, 5.14] 1.53* [ 1.03, 2.02] 0.63* [ 0.34, 0.92]
is_B 25.57* [18.34, 35.70] 101.65* [ 89.08, 116.52] 4.56* [ 4.32, 4.81] 1.42* [ 0.91, 1.92] 0.44* [ 0.18, 0.72]
Within-Person Effects
is_A:persuasion_self_cw 1.25* [ 1.07, 1.48] 1.00 [ 0.97, 1.04] -0.01 [-0.05, 0.04] 0.07 [-0.03, 0.18] -0.07 [-0.19, 0.04]
is_B:persuasion_self_cw 1.25* [ 1.08, 1.46] 1.05* [ 1.01, 1.09] 0.01 [-0.05, 0.06] 0.04 [-0.07, 0.16] 0.00 [-0.08, 0.09]
is_A:persuasion_partner_cw 1.24* [ 1.07, 1.45] 1.03 [ 0.99, 1.07] 0.01 [-0.04, 0.06] 0.07 [-0.02, 0.16] 0.09 [ 0.00, 0.18]
is_B:persuasion_partner_cw 1.28* [ 1.10, 1.55] 1.02 [ 0.98, 1.05] 0.05 [-0.01, 0.10] 0.08 [-0.01, 0.16] -0.06 [-0.15, 0.04]
is_A:day 0.66 [ 0.43, 1.04] 0.94 [ 0.83, 1.06] 0.08 [-0.10, 0.26] -0.01 [-0.39, 0.39] -0.04 [-0.42, 0.31]
is_B:day 0.87 [ 0.55, 1.37] 1.00 [ 0.89, 1.12] 0.40* [ 0.20, 0.58] -0.35 [-0.74, 0.04] 0.32 [-0.08, 0.71]
is_A:weartime_self_cw NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cw NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb 0.88 [ 0.43, 1.85] 0.88 [ 0.68, 1.13] 0.54 [-0.15, 1.22] -0.38 [-1.53, 0.83] 0.51 [-0.06, 1.10]
is_B:persuasion_self_cb 1.45 [ 0.72, 2.95] 1.28 [ 0.90, 1.79] -0.50 [-1.22, 0.20] -0.41 [-1.96, 0.98] 0.56* [ 0.05, 1.05]
is_A:persuasion_partner_cb 1.14 [ 0.50, 2.66] 1.29 [ 0.97, 1.72] -0.40 [-1.24, 0.45] 1.16 [-0.16, 2.55] 0.37 [-0.26, 1.01]
is_B:persuasion_partner_cb 0.60 [ 0.32, 1.14] 0.88 [ 0.66, 1.18] 0.51 [-0.08, 1.11] 0.92 [-0.30, 2.27] 0.29 [-0.14, 0.73]
is_A:weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(is_A) 0.78 [ 0.56, 1.05] 0.28 [0.21, 0.37] 0.79 [0.60, 1.03] 1.40 [1.09, 1.83] 0.59 [ 0.38, 0.86]
sd(is_B) 0.63 [ 0.44, 0.86] 0.34 [0.26, 0.44] 0.70 [0.54, 0.91] 1.46 [1.14, 1.87] 0.41 [ 0.21, 0.66]
sd(is_A:persuasion_self_cw) 0.15 [ 0.01, 0.33] 0.06 [0.01, 0.10] 0.05 [0.00, 0.13] 0.14 [0.01, 0.28] 0.23 [ 0.11, 0.37]
sd(is_B:persuasion_self_cw) 0.11 [ 0.01, 0.30] 0.06 [0.02, 0.11] 0.08 [0.02, 0.16] 0.21 [0.08, 0.35] 0.09 [ 0.01, 0.21]
sd(is_A:persuasion_partner_cw) 0.08 [ 0.00, 0.24] 0.06 [0.01, 0.12] 0.06 [0.00, 0.13] 0.07 [0.00, 0.19] 0.09 [ 0.01, 0.20]
sd(is_B:persuasion_partner_cw) 0.18 [ 0.02, 0.39] 0.05 [0.01, 0.10] 0.10 [0.03, 0.17] 0.05 [0.00, 0.16] 0.09 [ 0.01, 0.21]
Additional Parameters
ar[1] 0.03 [-0.93, 0.93] 0.23 [0.19, 0.26] 0.33 [0.29, 0.36] 0.33 [0.30, 0.37] -0.03 [-0.12, 0.06]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.14 [ 0.13, 0.15] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.05 [ 0.00, 0.13] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.54 [0.53, 0.55] 0.83 [0.81, 0.85] 1.74 [1.70, 1.78] 0.95 [ 0.90, 1.01]
# prepare export to excel
all_models$Persuasion_apim <- prepare_for_export(persuasion_all_apim_models)

Modelling Pressure

# For indistinguishable Dyads
model_rows_fixed_pressure_ind <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'day', 
    'weartime_self_cw', 
    'barriers_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'pressure_self_cb',
    'pressure_partner_cb',
    'weartime_self_cb',
    'barriers_self_cb'
  )
  
model_rows_random_pressure_ind <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)


# For APIMs

model_rows_fixed_pressure_apim <- c(
    'is_A',
    'is_B',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'is_A:pressure_self_cw',
    'is_B:pressure_self_cw', 
    'is_A:pressure_partner_cw', 
    'is_B:pressure_partner_cw',
    'is_A:day', 
    'is_B:day', 
    'is_A:weartime_self_cw', 
    'is_B:weartime_self_cw', 
    'is_A:barriers_self_cw',
    'is_B:barriers_self_cw',
    
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'is_A:pressure_self_cb',
    'is_B:pressure_self_cb', 
    'is_A:pressure_partner_cb', 
    'is_B:pressure_partner_cb',
    'is_A:weartime_self_cb', 
    'is_B:weartime_self_cb', 
    'is_A:barriers_self_cb',
    'is_B:barriers_self_cb'
  )
  

model_rows_random_pressure_apim <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(is_A)',
  'sd(is_B)',
  'sd(is_A:pressure_self_cw)',
  'sd(is_B:pressure_self_cw)', 
  'sd(is_A:pressure_partner_cw)', 
  'sd(is_B:pressure_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

Subjective MVPA ~ pressure

Indistinguishable

formula <- bf(
  pa_sub ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day + 
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_sub_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pressure_pa_sub_ind"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(pressure_pa_sub_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pressure_pa_sub_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_pa_sub_ind)
## Warning: Found 1 observations with a pareto_k > 0.7 in
## model 'pressure_pa_sub_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
print(loo_ind)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12080.2 177.4
## p_loo        20.3   1.9
## looic     24160.5 354.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 2.7]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3735  100.0%  260     
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_pa_sub_ind, ask = FALSE)

summarize_brms(
  pressure_pa_sub_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 32.32* 24.51 42.45 1.000 2554.90 2982.96
Within-Person Effects
pressure_self_cw 1.10 0.85 1.46 1.001 6708.20 3364.90
pressure_partner_cw 1.29 0.99 1.77 1.000 6458.10 2655.47
day 0.62* 0.45 0.85 1.001 8969.83 2796.94
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb 0.97 0.61 1.56 1.002 6241.32 2769.63
pressure_partner_cb 0.71 0.44 1.13 1.002 6036.00 2712.48
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.63 0.47 0.84 1.00 1490.47 2249.58
sd(pressure_self_cw) 0.16 0.01 0.48 1.00 2328.92 2325.71
sd(pressure_partner_cw) 0.16 0.00 0.49 1.00 2444.19 1688.36
Additional Parameters
ar[1] 0.03 -0.93 0.94 1.00 6953.75 2741.97
nu NA NA NA NA NA NA
shape 0.13 0.13 0.14 1.00 9762.02 2834.95
sderr 0.05 0.00 0.13 1.00 3080.77 2550.99
sigma NA NA NA NA NA NA

APIM

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_sub_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pressure_pa_sub_apim"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(pressure_pa_sub_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pressure_pa_sub_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_pa_sub_apim)
## Warning: Found 6 observations with a pareto_k > 0.7 in
## model 'pressure_pa_sub_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12079.6 177.5
## p_loo        31.6   2.7
## looic     24159.3 354.9
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.3]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3730  99.8%   519     
##    (0.7, 1]   (bad)         6   0.2%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_pa_sub_apim, ask = FALSE)

summarize_brms(
  pressure_pa_sub_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 36.54* 26.06 51.27 1.002 2413.33 3056.26
is_B 29.66* 21.84 41.14 1.002 2886.16 2801.38
Within-Person Effects
is_A:pressure_self_cw 1.05 0.72 1.59 1.002 4466.04 2441.66
is_B:pressure_self_cw 1.41 0.81 3.16 1.002 3498.53 2613.40
is_A:pressure_partner_cw 1.34 0.88 2.23 1.000 4367.48 3113.23
is_B:pressure_partner_cw 1.48 0.96 2.59 1.001 5029.86 3182.85
is_A:day 0.50* 0.33 0.77 1.001 7097.40 2963.09
is_B:day 0.69 0.45 1.07 1.001 4753.90 2773.37
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb 0.77 0.24 2.41 1.000 1986.14 2393.89
is_B:pressure_self_cb 1.04 0.29 3.89 1.000 2224.84 2834.00
is_A:pressure_partner_cb 1.06 0.24 4.75 1.001 2043.94 2311.01
is_B:pressure_partner_cb 0.51 0.17 1.46 1.000 2050.53 2482.69
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.74 0.53 1.00 1.00 2057.74 2889.82
sd(is_B) 0.60 0.41 0.83 1.00 2200.07 2983.87
sd(is_A:pressure_self_cw) 0.28 0.01 0.97 1.00 1841.59 2480.61
sd(is_B:pressure_self_cw) 0.47 0.02 1.52 1.00 1803.49 2223.63
sd(is_A:pressure_partner_cw) 0.28 0.01 0.94 1.00 2794.58 2621.34
sd(is_B:pressure_partner_cw) 0.32 0.01 0.94 1.00 2313.95 2563.18
Additional Parameters
ar[1] 0.02 -0.92 0.93 1.00 6166.03 2692.45
nu NA NA NA NA NA NA
shape 0.13 0.13 0.14 1.00 8383.44 2728.27
sderr 0.05 0.00 0.13 1.00 2713.49 2135.62
sigma NA NA NA NA NA NA

Compare models

loo_compare(loo_ind, loo_apim)
##                      elpd_diff se_diff
## pressure_pa_sub_apim  0.0       0.0   
## pressure_pa_sub_ind  -0.6       2.8
bayes_R2(pressure_pa_sub_apim)
##     Estimate Est.Error      Q2.5    Q97.5
## R2 0.3309817 0.1320252 0.1238244 0.500204
bayes_R2(pressure_pa_sub_ind)
##     Estimate Est.Error       Q2.5     Q97.5
## R2 0.1561202 0.0688205 0.08302701 0.3799566

Device Based MVPA ~ pressure

Indistinguishable

formula <- bf(
  pa_obj_log ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_obj_log_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pressure_pa_obj_log_ind"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(pressure_pa_obj_log_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pressure_pa_obj_log_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_pa_obj_log_ind)
print(loo_ind)
## 
## Computed from 4000 by 3299 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2815.0  55.5
## p_loo        51.5   2.4
## looic      5629.9 111.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 2.3]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pressure_pa_obj_log_ind, ask = FALSE)

summarize_brms(
  pressure_pa_obj_log_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
## Warning in summarize_brms(pressure_pa_obj_log_ind,
## model_rows_fixed = model_rows_fixed_pressure_ind, :
## Coefficients were exponentiated. Double check if this
## was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 120.04* 108.21 132.85 1.003 668.56 1463.19
Within-Person Effects
pressure_self_cw 0.99 0.92 1.05 1.002 4035.76 2618.01
pressure_partner_cw 1.02 0.96 1.08 1.001 5384.18 2778.75
day 0.93 0.85 1.02 1.000 6526.72 2723.65
weartime_self_cw 1.00* 1.00 1.00 1.001 4150.83 2457.88
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb 0.98 0.81 1.17 1.001 3697.25 2919.80
pressure_partner_cb 1.08 0.91 1.26 1.000 3464.73 2795.87
weartime_self_cb 1.00 1.00 1.00 1.001 3250.61 3145.83
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.29 0.22 0.38 1.01 821.24 1537.26
sd(pressure_self_cw) 0.05 0.00 0.15 1.00 1467.31 1843.60
sd(pressure_partner_cw) 0.04 0.00 0.11 1.00 2109.95 1766.72
Additional Parameters
ar[1] 0.29 0.26 0.32 1.00 5744.31 3205.59
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.56 0.55 0.58 1.00 7482.83 2899.76

apim

formula <- bf(
  pa_obj_log ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_obj_log_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pressure_pa_obj_log_apim"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(pressure_pa_obj_log_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pressure_pa_obj_log_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_pa_obj_log_apim)
## Warning: Found 4 observations with a pareto_k > 0.7 in
## model 'pressure_pa_obj_log_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 3299 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2742.8  56.7
## p_loo        95.2   5.5
## looic      5485.7 113.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3295  99.9%   321     
##    (0.7, 1]   (bad)         4   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_pa_obj_log_apim, ask = FALSE)

summarize_brms(
  pressure_pa_obj_log_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
## Warning in summarize_brms(pressure_pa_obj_log_apim,
## model_rows_fixed = model_rows_fixed_pressure_apim, :
## Coefficients were exponentiated. Double check if this
## was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 138.47* 123.18 155.40 1.001 1213.96 2060.96
is_B 102.06* 89.80 116.38 1.001 1339.98 2268.92
Within-Person Effects
is_A:pressure_self_cw 0.99 0.91 1.07 1.000 5105.44 2873.19
is_B:pressure_self_cw 0.98 0.86 1.10 1.000 2807.78 2656.54
is_A:pressure_partner_cw 1.07 0.92 1.27 1.001 3075.40 3028.11
is_B:pressure_partner_cw 1.00 0.93 1.09 1.001 5855.75 2878.22
is_A:day 0.91 0.81 1.02 1.002 5513.42 3009.55
is_B:day 0.96 0.86 1.08 1.000 5561.03 3014.30
is_A:weartime_self_cw 1.00* 1.00 1.00 1.003 3601.44 2581.10
is_B:weartime_self_cw 1.00* 1.00 1.00 1.001 4045.01 2852.39
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb 0.90 0.57 1.47 1.002 1363.63 1685.80
is_B:pressure_self_cb 0.88 0.44 1.68 1.006 1005.39 1443.61
is_A:pressure_partner_cb 1.21 0.67 2.22 1.003 1304.34 1764.46
is_B:pressure_partner_cb 1.19 0.71 2.03 1.005 1144.61 1935.66
is_A:weartime_self_cb 1.00 1.00 1.00 1.001 1286.41 2140.78
is_B:weartime_self_cb 1.00 1.00 1.00 1.002 2217.06 2726.55
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.29 0.22 0.38 1.00 1206.71 2206.69
sd(is_B) 0.34 0.26 0.45 1.00 1326.50 2330.11
sd(is_A:pressure_self_cw) 0.05 0.00 0.15 1.00 2159.15 2259.81
sd(is_B:pressure_self_cw) 0.11 0.01 0.28 1.00 1402.99 1819.43
sd(is_A:pressure_partner_cw) 0.19 0.01 0.42 1.00 1203.18 1228.82
sd(is_B:pressure_partner_cw) 0.04 0.00 0.13 1.00 2540.88 1920.82
Additional Parameters
ar[1] 0.22 0.19 0.26 1.00 5955.10 3249.08
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.55 0.53 0.56 1.00 7496.81 2898.87

Compare Models

loo_compare(loo_ind, loo_apim)
##                          elpd_diff se_diff
## pressure_pa_obj_log_apim   0.0       0.0  
## pressure_pa_obj_log_ind  -72.1      14.0
bayes_R2(pressure_pa_obj_log_apim)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.1834497 0.006406561 0.1707924 0.1959675
bayes_R2(pressure_pa_obj_log_ind)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.1633216 0.00669006 0.1502936 0.1766312

Affect ~ pressure

Indistinguishable

formula <- bf(
  aff ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_aff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pressure_aff_ind"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(pressure_aff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pressure_aff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_aff_ind)
print(loo_ind)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4821.5  63.4
## p_loo        58.9   3.6
## looic      9642.9 126.9
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.5]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pressure_aff_ind, ask = FALSE)

summarize_brms(
  pressure_aff_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.74* 4.52 4.96 1.004 851.16 1451.11
Within-Person Effects
pressure_self_cw -0.03 -0.14 0.06 1.002 3604.88 2917.81
pressure_partner_cw 0.02 -0.10 0.13 1.001 3155.97 2241.77
day 0.20* 0.04 0.36 1.001 11656.82 3091.57
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb 0.04 -0.26 0.34 1.000 5404.62 3390.21
pressure_partner_cb 0.04 -0.27 0.35 1.000 6114.32 2938.14
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.61 0.48 0.79 1.00 959.39 1913.84
sd(pressure_self_cw) 0.11 0.01 0.29 1.01 1276.44 1849.10
sd(pressure_partner_cw) 0.16 0.01 0.35 1.00 1037.86 1445.63
Additional Parameters
ar[1] 0.45 0.42 0.48 1.00 7934.28 3016.94
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.87 0.85 0.89 1.00 9434.80 2983.53

APIM

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_aff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pressure_aff_apim"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(pressure_aff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pressure_aff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_aff_apim)
## Warning: Found 3 observations with a pareto_k > 0.7 in
## model 'pressure_aff_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4684.1  65.4
## p_loo        99.8   4.7
## looic      9368.1 130.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.6]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3733  99.9%   99      
##    (0.7, 1]   (bad)         3   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_aff_apim, ask = FALSE)

summarize_brms(
  pressure_aff_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 4.85* 4.56 5.14 1.014 422.93 902.84
is_B 4.58* 4.31 4.84 1.004 745.99 1471.26
Within-Person Effects
is_A:pressure_self_cw -0.01 -0.13 0.11 1.000 4396.27 2881.27
is_B:pressure_self_cw -0.08 -0.23 0.06 1.000 3854.43 2376.21
is_A:pressure_partner_cw -0.03 -0.25 0.18 1.002 2974.97 2995.08
is_B:pressure_partner_cw 0.05 -0.07 0.17 1.003 5527.66 2827.12
is_A:day 0.07 -0.11 0.25 1.001 6238.35 2813.61
is_B:day 0.36* 0.18 0.55 1.002 7903.77 2504.67
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb 0.25 -0.97 1.40 1.001 1417.88 1791.06
is_B:pressure_self_cb -0.28 -1.69 1.08 1.001 1528.57 2066.28
is_A:pressure_partner_cb -0.23 -1.72 1.29 1.001 1444.00 2030.80
is_B:pressure_partner_cb 0.30 -0.75 1.41 1.001 1560.32 2088.62
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.81 0.63 1.05 1.00 845.83 1526.39
sd(is_B) 0.71 0.55 0.94 1.00 964.33 1499.46
sd(is_A:pressure_self_cw) 0.12 0.00 0.34 1.00 1031.55 1728.95
sd(is_B:pressure_self_cw) 0.14 0.00 0.40 1.00 1546.95 2185.63
sd(is_A:pressure_partner_cw) 0.29 0.03 0.67 1.00 1136.14 1168.30
sd(is_B:pressure_partner_cw) 0.09 0.00 0.27 1.00 2053.93 1960.06
Additional Parameters
ar[1] 0.33 0.30 0.36 1.00 7139.36 2988.17
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.84 0.82 0.86 1.00 7719.09 2923.68

Compare Models

loo_compare(loo_ind, loo_apim)
##                   elpd_diff se_diff
## pressure_aff_apim    0.0       0.0 
## pressure_aff_ind  -137.4      19.4
bayes_R2(pressure_aff_apim)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2363495 0.004612259 0.2271535 0.2452849
bayes_R2(pressure_aff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2257496 0.004834356 0.2161791 0.2349657

Anticipatory Affect ~ pressure

Indistinguishable

formula <- bf(
  anticipAff ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=-5, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_anticipAff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pressure_anticipAff_ind"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(pressure_anticipAff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pressure_anticipAff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_anticipAff_ind)
print(loo_ind)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7519.8  56.3
## p_loo        53.2   2.3
## looic     15039.5 112.6
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 2.0]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pressure_anticipAff_ind, ask = FALSE)

summarize_brms(
  pressure_anticipAff_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.50* 1.03 1.97 1.004 643.07 896.49
Within-Person Effects
pressure_self_cw 0.14 -0.10 0.34 1.001 2963.12 2621.66
pressure_partner_cw 0.28* 0.09 0.45 1.002 3020.37 2350.68
day -0.24 -0.57 0.09 1.001 6078.36 2705.64
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb -0.06 -0.67 0.53 1.001 3500.78 2425.01
pressure_partner_cb 0.51 -0.09 1.13 1.000 3318.45 2756.81
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.26 0.99 1.62 1.00 735.76 1305.28
sd(pressure_self_cw) 0.24 0.03 0.51 1.00 1417.58 1503.62
sd(pressure_partner_cw) 0.13 0.01 0.37 1.00 1514.30 1895.10
Additional Parameters
ar[1] 0.42 0.39 0.46 1.00 5599.05 2830.55
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.80 1.76 1.84 1.00 5283.11 2851.80

APIM

formula <- bf(
  anticipAff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_anticipAff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pressure_anticipAff_apim"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(pressure_anticipAff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pressure_anticipAff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_anticipAff_apim)
print(loo_apim)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7417.8  56.2
## p_loo        95.6   3.9
## looic     14835.6 112.5
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.4]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pressure_anticipAff_apim, ask = FALSE)

summarize_brms(
  pressure_anticipAff_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 1.55* 1.04 2.07 1.004 738.46 1503.76
is_B 1.47* 0.91 2.01 1.001 883.16 1596.89
Within-Person Effects
is_A:pressure_self_cw 0.19 -0.16 0.48 1.000 3360.55 2566.04
is_B:pressure_self_cw 0.04 -0.35 0.37 1.000 3316.27 2952.54
is_A:pressure_partner_cw 0.18 -0.18 0.50 1.001 2673.29 3399.08
is_B:pressure_partner_cw 0.27* 0.02 0.52 0.999 4412.25 2763.79
is_A:day -0.08 -0.47 0.31 1.001 7531.58 2493.13
is_B:day -0.39* -0.77 0.00 1.001 6513.32 3136.49
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb -0.51 -2.61 1.69 1.002 1335.12 2010.76
is_B:pressure_self_cb 0.32 -2.50 3.14 1.005 1358.05 1973.70
is_A:pressure_partner_cb 1.12 -1.69 3.76 1.003 1324.68 1860.98
is_B:pressure_partner_cb 0.27 -2.04 2.56 1.004 1344.82 2037.20
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 1.45 1.14 1.85 1.00 1121.97 1812.02
sd(is_B) 1.50 1.17 1.93 1.00 1194.02 2286.66
sd(is_A:pressure_self_cw) 0.32 0.02 0.73 1.01 1146.93 1060.96
sd(is_B:pressure_self_cw) 0.31 0.02 0.70 1.00 1466.16 1031.91
sd(is_A:pressure_partner_cw) 0.32 0.02 0.75 1.00 1833.35 1955.23
sd(is_B:pressure_partner_cw) 0.17 0.01 0.49 1.00 2000.75 2495.65
Additional Parameters
ar[1] 0.34 0.30 0.37 1.00 7363.28 3124.59
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.74 1.70 1.78 1.00 9870.15 2767.09

Compare Models

loo_compare(loo_ind, loo_apim)
##                          elpd_diff se_diff
## pressure_anticipAff_apim    0.0       0.0 
## pressure_anticipAff_ind  -101.9      15.8
bayes_R2(pressure_anticipAff_apim)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2200751 0.00455847 0.2109987 0.2287366
bayes_R2(pressure_anticipAff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2114515 0.004861903 0.2021107 0.2210266

reactance ~ pressure

Indistinguishable

formula <- bf(
  reactance ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_reactance_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  control = list(adapt_delta = 0.90),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pressure_reactance_ind"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(pressure_reactance_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pressure_reactance_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_reactance_ind)
## Warning: Found 5 observations with a pareto_k > 0.7 in
## model 'pressure_reactance_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
print(loo_ind)
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1055.0 35.0
## p_loo        50.3  5.9
## looic      2110.0 70.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.7]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     751   99.3%   106     
##    (0.7, 1]   (bad)        5    0.7%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_reactance_ind, ask = FALSE)

summarize_brms(
  pressure_reactance_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
## Warning: There were 1 divergent transitions after
## warmup. Increasing adapt_delta above 0.9 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.45* 0.30 0.59 1.000 3193.74 3437.17
Within-Person Effects
pressure_self_cw 0.31* 0.11 0.50 1.003 1611.45 2370.22
pressure_partner_cw 0.14 -0.04 0.37 1.001 1731.70 2125.96
day 0.16 -0.08 0.41 1.001 4758.13 2572.79
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb 0.65* 0.41 0.89 1.000 3036.44 3152.13
pressure_partner_cb 0.11 -0.13 0.35 1.001 3160.95 3068.98
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.26 0.16 0.38 1.00 2049.13 2273.94
sd(pressure_self_cw) 0.36 0.20 0.58 1.00 1372.23 2312.18
sd(pressure_partner_cw) 0.21 0.01 0.56 1.00 767.61 1434.60
Additional Parameters
ar[1] 0.00 -0.08 0.08 1.00 4270.81 2948.97
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.94 0.89 0.99 1.00 3866.73 2892.88

APIM

formula <- bf(
  reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 

  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_reactance_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  control = list(adapt_delta = 0.90, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pressure_reactance_apim"
)
## Warning: Rows containing NAs were excluded from the
## model.
pp_check(pressure_reactance_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pressure_reactance_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_reactance_apim)
## Warning: Found 16 observations with a pareto_k > 0.7 in
## model 'pressure_reactance_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1052.6 35.2
## p_loo        68.8  6.5
## looic      2105.2 70.3
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     740   97.9%   34      
##    (0.7, 1]   (bad)       15    2.0%   <NA>    
##    (1, Inf)   (very bad)   1    0.1%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_reactance_apim, ask = FALSE)

summarize_brms(
  pressure_reactance_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 0.47* 0.27 0.66 1.001 3865.17 2625.40
is_B 0.43* 0.23 0.64 1.000 4639.50 3231.41
Within-Person Effects
is_A:pressure_self_cw 0.28* 0.04 0.52 1.000 2549.46 2660.78
is_B:pressure_self_cw 0.35* 0.06 0.67 1.001 2771.46 2812.13
is_A:pressure_partner_cw 0.24 -0.12 0.60 1.001 3444.68 3041.04
is_B:pressure_partner_cw 0.01 -0.17 0.21 1.001 4982.60 2800.85
is_A:day 0.02 -0.29 0.36 1.002 6111.56 3255.57
is_B:day 0.25 -0.13 0.60 1.000 4600.77 2778.98
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb 0.66* 0.03 1.29 1.002 2636.99 3015.18
is_B:pressure_self_cb 0.68 -0.02 1.43 1.000 2602.90 2540.28
is_A:pressure_partner_cb 0.05 -0.76 0.86 1.001 2771.70 3021.92
is_B:pressure_partner_cb 0.07 -0.50 0.63 1.001 2837.72 2685.37
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.33 0.19 0.49 1.00 1753.48 2973.59
sd(is_B) 0.31 0.17 0.50 1.00 2175.49 2977.66
sd(is_A:pressure_self_cw) 0.35 0.16 0.63 1.00 2131.44 2373.76
sd(is_B:pressure_self_cw) 0.45 0.18 0.83 1.00 1915.06 1887.59
sd(is_A:pressure_partner_cw) 0.34 0.02 0.84 1.00 1579.62 1857.61
sd(is_B:pressure_partner_cw) 0.12 0.00 0.45 1.00 2117.20 2364.57
Additional Parameters
ar[1] -0.03 -0.12 0.05 1.00 5620.36 3236.28
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.92 0.88 0.97 1.00 5363.64 3014.61

Compare Models

loo_compare(loo_ind, loo_apim)
##                         elpd_diff se_diff
## pressure_reactance_apim  0.0       0.0   
## pressure_reactance_ind  -2.4       4.8
bayes_R2(pressure_reactance_apim)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2417501 0.009366639 0.2232442 0.2593775
bayes_R2(pressure_reactance_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2344929 0.009454708 0.2150926 0.2527116

Report models for pressure side by side

Indistinguishable Models

pressure_all_ind_models <- report_side_by_side(
  pressure_pa_sub_ind,
  pressure_pa_obj_log_ind,
  pressure_aff_ind,
  pressure_anticipAff_ind,
  pressure_reactance_ind,
  
  model_rows_random = model_rows_random_pressure_ind,
  model_rows_fixed = model_rows_fixed_pressure_ind
) 
## [1] "pressure_pa_sub_ind"
## [1] "pressure_pa_obj_log_ind"
## Warning in summarize_brms(model, short_version = TRUE,
## exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "pressure_aff_ind"
## [1] "pressure_anticipAff_ind"
## [1] "pressure_reactance_ind"
## Warning: There were 1 divergent transitions after
## warmup. Increasing adapt_delta above 0.9 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
# pretty printing
pressure_all_ind_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_ind()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR pressure_pa_sub_ind 95% CI pressure_pa_sub_ind exp(Est.) pressure_pa_obj_log_ind 95% CI pressure_pa_obj_log_ind b pressure_aff_ind 95% CI pressure_aff_ind b pressure_anticipAff_ind 95% CI pressure_anticipAff_ind b pressure_reactance_ind 95% CI pressure_reactance_ind
Intercept 32.32* [24.51, 42.45] 120.04* [108.21, 132.85] 4.74* [ 4.52, 4.96] 1.50* [ 1.03, 1.97] 0.45* [ 0.30, 0.59]
Within-Person Effects
pressure_self_cw 1.10 [ 0.85, 1.46] 0.99 [ 0.92, 1.05] -0.03 [-0.14, 0.06] 0.14 [-0.10, 0.34] 0.31* [ 0.11, 0.50]
pressure_partner_cw 1.29 [ 0.99, 1.77] 1.02 [ 0.96, 1.08] 0.02 [-0.10, 0.13] 0.28* [ 0.09, 0.45] 0.14 [-0.04, 0.37]
day 0.62* [ 0.45, 0.85] 0.93 [ 0.85, 1.02] 0.20* [ 0.04, 0.36] -0.24 [-0.57, 0.09] 0.16 [-0.08, 0.41]
weartime_self_cw NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb 0.97 [ 0.61, 1.56] 0.98 [ 0.81, 1.17] 0.04 [-0.26, 0.34] -0.06 [-0.67, 0.53] 0.65* [ 0.41, 0.89]
pressure_partner_cb 0.71 [ 0.44, 1.13] 1.08 [ 0.91, 1.26] 0.04 [-0.27, 0.35] 0.51 [-0.09, 1.13] 0.11 [-0.13, 0.35]
weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.63 [ 0.47, 0.84] 0.29 [0.22, 0.38] 0.61 [0.48, 0.79] 1.26 [0.99, 1.62] 0.26 [ 0.16, 0.38]
sd(pressure_self_cw) 0.16 [ 0.01, 0.48] 0.05 [0.00, 0.15] 0.11 [0.01, 0.29] 0.24 [0.03, 0.51] 0.36 [ 0.20, 0.58]
sd(pressure_partner_cw) 0.16 [ 0.00, 0.49] 0.04 [0.00, 0.11] 0.16 [0.01, 0.35] 0.13 [0.01, 0.37] 0.21 [ 0.01, 0.56]
Additional Parameters
ar[1] 0.03 [-0.93, 0.94] 0.29 [0.26, 0.32] 0.45 [0.42, 0.48] 0.42 [0.39, 0.46] 0.00 [-0.08, 0.08]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.13 [ 0.13, 0.14] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.05 [ 0.00, 0.13] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.56 [0.55, 0.58] 0.87 [0.85, 0.89] 1.80 [1.76, 1.84] 0.94 [ 0.89, 0.99]
# prepare export to excel
all_models$Pressure_ind <- prepare_for_export(pressure_all_ind_models)

APIMs

pressure_all_apim_models <- report_side_by_side(
  pressure_pa_sub_apim,
  pressure_pa_obj_log_apim,
  pressure_aff_apim,
  pressure_anticipAff_apim,
  pressure_reactance_apim,
  
  model_rows_random = model_rows_random_pressure_apim,
  model_rows_fixed = model_rows_fixed_pressure_apim
) 
## [1] "pressure_pa_sub_apim"
## [1] "pressure_pa_obj_log_apim"
## Warning in summarize_brms(model, short_version = TRUE,
## exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "pressure_aff_apim"
## [1] "pressure_anticipAff_apim"
## [1] "pressure_reactance_apim"
# pretty printing
pressure_all_apim_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_apim()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR pressure_pa_sub_apim 95% CI pressure_pa_sub_apim exp(Est.) pressure_pa_obj_log_apim 95% CI pressure_pa_obj_log_apim b pressure_aff_apim 95% CI pressure_aff_apim b pressure_anticipAff_apim 95% CI pressure_anticipAff_apim b pressure_reactance_apim 95% CI pressure_reactance_apim
is_A 36.54* [26.06, 51.27] 138.47* [123.18, 155.40] 4.85* [ 4.56, 5.14] 1.55* [ 1.04, 2.07] 0.47* [ 0.27, 0.66]
is_B 29.66* [21.84, 41.14] 102.06* [ 89.80, 116.38] 4.58* [ 4.31, 4.84] 1.47* [ 0.91, 2.01] 0.43* [ 0.23, 0.64]
Within-Person Effects
is_A:pressure_self_cw 1.05 [ 0.72, 1.59] 0.99 [ 0.91, 1.07] -0.01 [-0.13, 0.11] 0.19 [-0.16, 0.48] 0.28* [ 0.04, 0.52]
is_B:pressure_self_cw 1.41 [ 0.81, 3.16] 0.98 [ 0.86, 1.10] -0.08 [-0.23, 0.06] 0.04 [-0.35, 0.37] 0.35* [ 0.06, 0.67]
is_A:pressure_partner_cw 1.34 [ 0.88, 2.23] 1.07 [ 0.92, 1.27] -0.03 [-0.25, 0.18] 0.18 [-0.18, 0.50] 0.24 [-0.12, 0.60]
is_B:pressure_partner_cw 1.48 [ 0.96, 2.59] 1.00 [ 0.93, 1.09] 0.05 [-0.07, 0.17] 0.27* [ 0.02, 0.52] 0.01 [-0.17, 0.21]
is_A:day 0.50* [ 0.33, 0.77] 0.91 [ 0.81, 1.02] 0.07 [-0.11, 0.25] -0.08 [-0.47, 0.31] 0.02 [-0.29, 0.36]
is_B:day 0.69 [ 0.45, 1.07] 0.96 [ 0.86, 1.08] 0.36* [ 0.18, 0.55] -0.39* [-0.77, 0.00] 0.25 [-0.13, 0.60]
is_A:weartime_self_cw NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cw NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb 0.77 [ 0.24, 2.41] 0.90 [ 0.57, 1.47] 0.25 [-0.97, 1.40] -0.51 [-2.61, 1.69] 0.66* [ 0.03, 1.29]
is_B:pressure_self_cb 1.04 [ 0.29, 3.89] 0.88 [ 0.44, 1.68] -0.28 [-1.69, 1.08] 0.32 [-2.50, 3.14] 0.68 [-0.02, 1.43]
is_A:pressure_partner_cb 1.06 [ 0.24, 4.75] 1.21 [ 0.67, 2.22] -0.23 [-1.72, 1.29] 1.12 [-1.69, 3.76] 0.05 [-0.76, 0.86]
is_B:pressure_partner_cb 0.51 [ 0.17, 1.46] 1.19 [ 0.71, 2.03] 0.30 [-0.75, 1.41] 0.27 [-2.04, 2.56] 0.07 [-0.50, 0.63]
is_A:weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(is_A) 0.74 [ 0.53, 1.00] 0.29 [0.22, 0.38] 0.81 [0.63, 1.05] 1.45 [1.14, 1.85] 0.33 [ 0.19, 0.49]
sd(is_B) 0.60 [ 0.41, 0.83] 0.34 [0.26, 0.45] 0.71 [0.55, 0.94] 1.50 [1.17, 1.93] 0.31 [ 0.17, 0.50]
sd(is_A:pressure_self_cw) 0.28 [ 0.01, 0.97] 0.05 [0.00, 0.15] 0.12 [0.00, 0.34] 0.32 [0.02, 0.73] 0.35 [ 0.16, 0.63]
sd(is_B:pressure_self_cw) 0.47 [ 0.02, 1.52] 0.11 [0.01, 0.28] 0.14 [0.00, 0.40] 0.31 [0.02, 0.70] 0.45 [ 0.18, 0.83]
sd(is_A:pressure_partner_cw) 0.28 [ 0.01, 0.94] 0.19 [0.01, 0.42] 0.29 [0.03, 0.67] 0.32 [0.02, 0.75] 0.34 [ 0.02, 0.84]
sd(is_B:pressure_partner_cw) 0.32 [ 0.01, 0.94] 0.04 [0.00, 0.13] 0.09 [0.00, 0.27] 0.17 [0.01, 0.49] 0.12 [ 0.00, 0.45]
Additional Parameters
ar[1] 0.02 [-0.92, 0.93] 0.22 [0.19, 0.26] 0.33 [0.30, 0.36] 0.34 [0.30, 0.37] -0.03 [-0.12, 0.05]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.13 [ 0.13, 0.14] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.05 [ 0.00, 0.13] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.55 [0.53, 0.56] 0.84 [0.82, 0.86] 1.74 [1.70, 1.78] 0.92 [ 0.88, 0.97]
# prepare export to excel
all_models$Pressure_apim <- prepare_for_export(pressure_all_apim_models)

Modelling Pushing

# For indistinguishable Dyads
model_rows_fixed_pushing_ind <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw', 
    'barriers_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    'barriers_self_cb'
  )
  
model_rows_random_pushing_ind <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)


# For APIMs

model_rows_fixed_pushing_apim <- c(
    'is_A',
    'is_B',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'is_A:pushing_self_cw',
    'is_B:pushing_self_cw', 
    'is_A:pushing_partner_cw', 
    'is_B:pushing_partner_cw',
    'is_A:day', 
    'is_B:day', 
    'is_A:weartime_self_cw', 
    'is_B:weartime_self_cw', 
    'is_A:barriers_self_cw',
    'is_B:barriers_self_cw',
    
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'is_A:pushing_self_cb',
    'is_B:pushing_self_cb', 
    'is_A:pushing_partner_cb', 
    'is_B:pushing_partner_cb',
    'is_A:weartime_self_cb', 
    'is_B:weartime_self_cb', 
    'is_A:barriers_self_cb',
    'is_B:barriers_self_cb'
  )
  

model_rows_random_pushing_apim <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(is_A)',
  'sd(is_B)',
  'sd(is_A:pushing_self_cw)',
  'sd(is_B:pushing_self_cw)', 
  'sd(is_A:pushing_partner_cw)', 
  'sd(is_B:pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

Subjective MVPA ~ pushing

Indistinguishable

formula <- bf(
  pa_sub ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day + 
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_sub_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pushing_pa_sub_ind"
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_pa_sub_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pushing_pa_sub_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_pa_sub_ind)
## Warning: Found 6 observations with a pareto_k > 0.7 in
## model 'pushing_pa_sub_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
print(loo_ind)
## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -6312.3  67.3
## p_loo        31.3   2.4
## looic     12624.6 134.6
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 2.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1334  99.6%   259     
##    (0.7, 1]   (bad)         6   0.4%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_pa_sub_ind, ask = FALSE)

summarize_brms(
  pushing_pa_sub_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 53.20* 39.52 70.18 1.001 1789.85 2188.78
Within-Person Effects
pushing_self_cw 1.05 0.88 1.21 1.004 2590.55 2294.26
pushing_partner_cw 1.06 0.91 1.23 1.004 2835.45 2598.84
day 0.66* 0.49 0.88 1.000 5781.03 2914.40
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw 1.10* 1.03 1.17 1.000 8325.96 3013.80
Between-Person Effects
pushing_self_cb 0.95 0.40 2.07 1.007 945.35 1670.65
pushing_partner_cb 0.72 0.30 1.60 1.006 907.49 1838.40
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb 1.12 0.96 1.31 1.001 3668.19 3314.36
Random Effects
sd(Intercept) 0.64 0.43 0.94 1.00 846.38 2007.38
sd(pushing_self_cw) 0.24 0.03 0.47 1.01 613.24 1061.85
sd(pushing_partner_cw) 0.21 0.02 0.45 1.02 461.51 855.82
Additional Parameters
ar[1] 0.04 -0.93 0.94 1.00 6135.24 2762.59
nu NA NA NA NA NA NA
shape 0.43 0.40 0.47 1.00 5970.29 2550.17
sderr 0.04 0.00 0.13 1.00 2616.00 2214.35
sigma NA NA NA NA NA NA

APIM

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_sub_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pushing_pa_sub_apim"
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_pa_sub_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pushing_pa_sub_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_pa_sub_apim)
## Warning: Found 6 observations with a pareto_k > 0.7 in
## model 'pushing_pa_sub_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -6317.1  67.4
## p_loo        43.8   2.9
## looic     12634.2 134.7
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.2, 1.5]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1334  99.6%   303     
##    (0.7, 1]   (bad)         6   0.4%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_pa_sub_apim, ask = FALSE)

summarize_brms(
  pushing_pa_sub_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 53.36* 36.08 78.05 1.002 1162.51 2005.31
is_B 51.56* 37.38 70.87 1.003 1721.93 2414.86
Within-Person Effects
is_A:pushing_self_cw 1.02 0.87 1.19 1.000 3124.23 2718.83
is_B:pushing_self_cw 1.09 0.81 1.38 1.005 1279.48 1211.04
is_A:pushing_partner_cw 1.11 0.87 1.38 1.002 1526.17 1290.88
is_B:pushing_partner_cw 1.06 0.91 1.23 1.000 4044.11 2857.58
is_A:day 0.58* 0.39 0.88 1.001 3526.13 3184.52
is_B:day 0.71 0.47 1.09 1.004 2597.27 2758.43
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw 1.06 0.97 1.17 1.000 4600.49 2909.64
is_B:barriers_self_cw 1.14* 1.05 1.24 1.000 4442.27 2678.55
Between-Person Effects
is_A:pushing_self_cb 0.85 0.18 3.81 1.004 1084.86 1960.81
is_B:pushing_self_cb 1.03 0.33 3.02 1.001 1710.99 2495.05
is_A:pushing_partner_cb 0.74 0.15 3.10 1.000 1524.32 2172.38
is_B:pushing_partner_cb 0.49 0.13 1.55 1.003 1200.31 1790.46
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb 1.13 0.84 1.52 1.002 1359.13 2042.66
is_B:barriers_self_cb 1.28 0.99 1.68 1.001 1746.64 2252.67
Random Effects
sd(is_A) 0.82 0.51 1.26 1.01 676.71 1724.27
sd(is_B) 0.55 0.34 0.82 1.00 1247.91 1906.65
sd(is_A:pushing_self_cw) 0.11 0.00 0.34 1.00 1604.98 1846.30
sd(is_B:pushing_self_cw) 0.34 0.03 0.80 1.01 478.90 1080.01
sd(is_A:pushing_partner_cw) 0.25 0.01 0.68 1.01 408.33 924.15
sd(is_B:pushing_partner_cw) 0.11 0.00 0.31 1.00 1396.07 1298.87
Additional Parameters
ar[1] 0.05 -0.93 0.95 1.00 3468.63 2387.90
nu NA NA NA NA NA NA
shape 0.44 0.40 0.48 1.00 3641.33 2509.80
sderr 0.05 0.00 0.13 1.00 1761.60 1475.04
sigma NA NA NA NA NA NA

Compare models

loo_compare(loo_ind, loo_apim)
##                     elpd_diff se_diff
## pushing_pa_sub_ind   0.0       0.0   
## pushing_pa_sub_apim -4.8       3.6
bayes_R2(pushing_pa_sub_apim)
##    Estimate  Est.Error      Q2.5     Q97.5
## R2 0.256768 0.06108899 0.1723171 0.4283762
bayes_R2(pushing_pa_sub_ind)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2165939 0.04688659 0.1467498 0.3265735

Device Based MVPA ~ pushing

Indistinguishable

formula <- bf(
  pa_obj_log ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day + weartime_self_cw + weartime_self_cb +
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_obj_log_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pushing_pa_obj_log_ind"
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_pa_obj_log_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pushing_pa_obj_log_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_pa_obj_log_ind)
## Warning: Found 1 observations with a pareto_k > 0.7 in
## model 'pushing_pa_obj_log_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
print(loo_ind)
## 
## Computed from 4000 by 1203 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -976.9 33.7
## p_loo        62.2  4.6
## looic      1953.9 67.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 2.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1202  99.9%   247     
##    (0.7, 1]   (bad)         1   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_pa_obj_log_ind, ask = FALSE)

summarize_brms(
  pushing_pa_obj_log_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
## Warning in summarize_brms(pushing_pa_obj_log_ind,
## model_rows_fixed = model_rows_fixed_pushing_ind, :
## Coefficients were exponentiated. Double check if this was
## intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 126.34* 110.67 143.56 1.004 1640.04 2783.68
Within-Person Effects
pushing_self_cw 1.03 0.97 1.09 1.001 2852.57 2586.01
pushing_partner_cw 1.03 0.99 1.07 1.000 5289.60 3240.89
day 0.92 0.81 1.05 1.003 7322.25 2825.69
weartime_self_cw 1.00* 1.00 1.00 1.001 4017.63 2437.17
barriers_self_cw 1.03* 1.01 1.05 1.002 8411.14 2387.77
Between-Person Effects
pushing_self_cb 0.99 0.70 1.39 1.001 2178.79 2753.60
pushing_partner_cb 1.10 0.76 1.57 1.004 2285.91 2635.72
weartime_self_cb 1.00 1.00 1.00 1.001 4716.99 3614.61
barriers_self_cb 1.04 0.97 1.12 1.002 4735.19 3210.98
Random Effects
sd(Intercept) 0.30 0.22 0.41 1.00 1614.45 2565.68
sd(pushing_self_cw) 0.10 0.02 0.18 1.01 616.98 462.68
sd(pushing_partner_cw) 0.04 0.00 0.10 1.00 1448.13 2532.33
Additional Parameters
ar[1] 0.27 0.21 0.34 1.00 4780.73 3110.93
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.53 0.51 0.55 1.00 4854.21 2824.58

apim

formula <- bf(
  pa_obj_log ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_obj_log_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pushing_pa_obj_log_apim"
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_pa_obj_log_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pushing_pa_obj_log_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_pa_obj_log_apim)
## Warning: Found 2 observations with a pareto_k > 0.7 in
## model 'pushing_pa_obj_log_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 1203 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -960.6 34.7
## p_loo        94.8  6.6
## looic      1921.2 69.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.6]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1201  99.8%   110     
##    (0.7, 1]   (bad)         2   0.2%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_pa_obj_log_apim, ask = FALSE)

summarize_brms(
  pushing_pa_obj_log_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
## Warning in summarize_brms(pushing_pa_obj_log_apim,
## model_rows_fixed = model_rows_fixed_pushing_apim, :
## Coefficients were exponentiated. Double check if this was
## intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 147.66* 126.44 171.48 1.002 1534.30 2134.13
is_B 108.34* 93.29 125.51 1.001 1705.49 2176.22
Within-Person Effects
is_A:pushing_self_cw 1.00 0.93 1.07 1.001 2761.66 2646.46
is_B:pushing_self_cw 1.07 0.99 1.16 1.000 2574.40 2785.81
is_A:pushing_partner_cw 1.07 1.00 1.15 1.001 3011.06 2701.25
is_B:pushing_partner_cw 0.99 0.94 1.05 1.001 4015.71 2977.91
is_A:day 0.89 0.75 1.07 1.001 3445.88 2762.23
is_B:day 0.98 0.82 1.18 1.000 3386.74 2844.89
is_A:weartime_self_cw 1.00 1.00 1.00 1.000 3919.56 2962.42
is_B:weartime_self_cw 1.00 1.00 1.00 1.000 3738.78 2578.10
is_A:barriers_self_cw 1.02 0.98 1.05 1.001 4298.30 3176.14
is_B:barriers_self_cw 1.04* 1.01 1.07 1.000 4858.03 3158.02
Between-Person Effects
is_A:pushing_self_cb 0.71 0.39 1.34 1.003 1376.78 1952.00
is_B:pushing_self_cb 1.10 0.61 2.01 1.001 1471.98 2121.84
is_A:pushing_partner_cb 1.43 0.80 2.56 1.001 1287.66 2042.83
is_B:pushing_partner_cb 0.95 0.50 1.79 1.003 1465.98 1901.63
is_A:weartime_self_cb 1.00 1.00 1.00 1.000 2167.83 2612.60
is_B:weartime_self_cb 1.00 1.00 1.00 1.000 2900.50 3165.74
is_A:barriers_self_cb 1.02 0.91 1.15 1.001 1499.54 1793.43
is_B:barriers_self_cb 1.12 0.97 1.28 1.001 2054.07 2477.98
Random Effects
sd(is_A) 0.31 0.22 0.43 1.00 1538.70 2439.93
sd(is_B) 0.31 0.22 0.43 1.00 1677.72 2550.23
sd(is_A:pushing_self_cw) 0.10 0.02 0.18 1.01 1087.99 1014.68
sd(is_B:pushing_self_cw) 0.11 0.02 0.21 1.00 905.60 517.49
sd(is_A:pushing_partner_cw) 0.08 0.00 0.17 1.00 951.01 805.38
sd(is_B:pushing_partner_cw) 0.03 0.00 0.09 1.00 2237.30 2314.64
Additional Parameters
ar[1] 0.22 0.15 0.28 1.00 3244.45 2954.41
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.51 0.49 0.54 1.00 3558.34 2789.14

Compare Models

loo_compare(loo_ind, loo_apim)
##                         elpd_diff se_diff
## pushing_pa_obj_log_apim   0.0       0.0  
## pushing_pa_obj_log_ind  -16.3       8.3
bayes_R2(pushing_pa_obj_log_apim)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2140487 0.01144776 0.1910864 0.2356184
bayes_R2(pushing_pa_obj_log_ind)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.1882489 0.01217251 0.1648978 0.2122497

Affect ~ pushing

Indistinguishable

formula <- bf(
  aff ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day +
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_aff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pushing_aff_ind"
)
## Warning: Rows containing NAs were excluded from the model.
## Warning: Chain 1 finished unexpectedly!
## Warning: Chain 3 finished unexpectedly!
## Warning: Chain 4 finished unexpectedly!
## Warning: 3 chain(s) finished unexpectedly!
## Warning: The returned fit object will only read in results
## of successful chains. Please use read_cmdstan_csv() to read
## the results of the failed chains separately.Use the
## $output(chain_id) method for more output of the failed
## chains.
pp_check(pushing_aff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pushing_aff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_aff_ind)
print(loo_ind)
## 
## Computed from 1000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1731.4 39.3
## p_loo        52.5  3.1
## looic      3462.7 78.5
## ------
## MCSE of elpd_loo is 0.3.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.4]).
## 
## All Pareto k estimates are good (k < 0.67).
## See help('pareto-k-diagnostic') for details.
plot(pushing_aff_ind, ask = FALSE)

summarize_brms(
  pushing_aff_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.82* 4.61 5.05 0.999 180.24 235.36
Within-Person Effects
pushing_self_cw 0.03 -0.03 0.10 1.001 577.49 594.77
pushing_partner_cw 0.09* 0.02 0.16 1.003 548.38 509.68
day 0.22* 0.02 0.42 1.001 1069.52 754.93
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw 0.12* 0.08 0.15 1.005 852.29 596.28
Between-Person Effects
pushing_self_cb 0.36 -0.28 0.96 1.000 219.35 275.01
pushing_partner_cb 0.23 -0.43 0.81 1.000 214.39 250.48
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb 0.27* 0.17 0.38 1.003 518.14 706.95
Random Effects
sd(Intercept) 0.61 0.47 0.78 1.01 224.47 476.26
sd(pushing_self_cw) 0.06 0.00 0.15 1.00 292.14 236.75
sd(pushing_partner_cw) 0.09 0.01 0.18 1.00 374.27 397.29
Additional Parameters
ar[1] 0.24 0.17 0.29 1.00 797.79 767.47
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.86 0.83 0.90 1.00 1252.18 821.88

APIM

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_aff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pushing_aff_apim"
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_aff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pushing_aff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_aff_apim)
print(loo_apim)
## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1701.8 39.9
## p_loo        90.9  5.2
## looic      3403.7 79.8
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.2]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pushing_aff_apim, ask = FALSE)

summarize_brms(
  pushing_aff_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 4.92* 4.62 5.21 1.003 1330.59 2092.19
is_B 4.68* 4.43 4.92 1.003 1870.02 2585.73
Within-Person Effects
is_A:pushing_self_cw -0.02 -0.10 0.07 1.001 4709.48 3021.14
is_B:pushing_self_cw 0.07 -0.04 0.17 1.002 3788.51 2370.12
is_A:pushing_partner_cw 0.12* 0.01 0.22 1.001 4072.38 2872.46
is_B:pushing_partner_cw 0.08 0.00 0.16 1.000 5393.09 3172.40
is_A:day 0.05 -0.20 0.30 1.000 6380.29 2517.59
is_B:day 0.38* 0.14 0.63 1.000 7084.66 2942.88
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw 0.12* 0.07 0.17 1.000 8805.35 2866.16
is_B:barriers_self_cw 0.13* 0.08 0.17 1.001 9132.52 3021.69
Between-Person Effects
is_A:pushing_self_cb 0.83 -0.54 2.17 1.001 1370.54 2088.62
is_B:pushing_self_cb -0.06 -1.03 0.95 1.001 2022.59 2509.56
is_A:pushing_partner_cb -0.17 -1.42 1.15 1.002 1739.06 2335.23
is_B:pushing_partner_cb 0.70 -0.39 1.80 1.001 1783.58 2304.43
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb 0.21 -0.08 0.49 1.003 1288.58 1802.34
is_B:barriers_self_cb 0.28 0.00 0.55 1.002 1711.94 2063.38
Random Effects
sd(is_A) 0.78 0.60 1.01 1.00 1160.59 2280.08
sd(is_B) 0.61 0.46 0.81 1.00 1483.62 2645.92
sd(is_A:pushing_self_cw) 0.08 0.00 0.21 1.00 1917.54 2365.56
sd(is_B:pushing_self_cw) 0.10 0.00 0.25 1.00 1162.55 1258.96
sd(is_A:pushing_partner_cw) 0.12 0.01 0.29 1.00 1316.74 1409.96
sd(is_B:pushing_partner_cw) 0.07 0.00 0.18 1.00 2019.03 2370.58
Additional Parameters
ar[1] 0.15 0.09 0.21 1.00 5369.03 2931.51
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.83 0.80 0.87 1.00 6707.97 2962.17

Compare Models

loo_compare(loo_ind, loo_apim)
##                  elpd_diff se_diff
## pushing_aff_apim   0.0       0.0  
## pushing_aff_ind  -29.5       9.9
bayes_R2(pushing_aff_apim)
##     Estimate   Est.Error      Q2.5   Q97.5
## R2 0.2140515 0.007457561 0.1994122 0.22818
bayes_R2(pushing_aff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.1992556 0.007741568 0.1842412 0.2140275

Anticipatory Affect ~ pushing

Indistinguishable

formula <- bf(
  anticipAff ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day + 
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=-5, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_anticipAff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pushing_anticipAff_ind"
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_anticipAff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pushing_anticipAff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_anticipAff_ind)
print(loo_ind)
## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -2750.9 29.5
## p_loo        49.9  2.4
## looic      5501.9 59.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.7]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pushing_anticipAff_ind, ask = FALSE)

summarize_brms(
  pushing_anticipAff_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.85* 1.39 2.30 1.002 1341.09 2084.56
Within-Person Effects
pushing_self_cw 0.00 -0.14 0.14 1.001 3442.05 2981.29
pushing_partner_cw 0.15* 0.01 0.31 1.001 3143.25 2484.79
day -0.25 -0.71 0.23 1.001 4825.75 3265.32
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw 0.18* 0.11 0.26 1.000 4610.94 2844.30
Between-Person Effects
pushing_self_cb 0.69 -0.55 1.90 1.000 1292.63 1893.60
pushing_partner_cb 0.87 -0.32 2.10 1.000 1305.14 1996.88
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb 0.56* 0.30 0.82 1.001 2520.15 2733.92
Random Effects
sd(Intercept) 1.09 0.82 1.43 1.00 1245.98 2068.35
sd(pushing_self_cw) 0.14 0.01 0.33 1.00 1303.04 1512.19
sd(pushing_partner_cw) 0.18 0.02 0.37 1.01 1057.27 1027.03
Additional Parameters
ar[1] 0.33 0.27 0.38 1.00 3638.04 2658.71
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.85 1.78 1.92 1.00 3545.11 3038.28

APIM

formula <- bf(
  anticipAff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_anticipAff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pushing_anticipAff_apim"
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_anticipAff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pushing_anticipAff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_anticipAff_apim)
print(loo_apim)
## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -2724.0 30.9
## p_loo        86.2  4.3
## looic      5448.0 61.8
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pushing_anticipAff_apim, ask = FALSE)

summarize_brms(
  pushing_anticipAff_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 1.88* 1.31 2.43 1.002 1672.00 2460.47
is_B 1.83* 1.24 2.39 1.003 1677.46 2700.54
Within-Person Effects
is_A:pushing_self_cw 0.07 -0.11 0.25 1.001 5955.15 2813.89
is_B:pushing_self_cw -0.04 -0.25 0.18 1.001 4804.71 2991.04
is_A:pushing_partner_cw 0.16 -0.06 0.39 1.001 4715.03 3347.93
is_B:pushing_partner_cw 0.13 -0.05 0.30 1.003 6259.68 2979.30
is_A:day -0.20 -0.78 0.37 1.001 7278.01 2882.16
is_B:day -0.31 -0.89 0.27 1.001 7525.11 3191.86
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw 0.10* 0.00 0.21 1.003 8858.58 3219.06
is_B:barriers_self_cw 0.26* 0.16 0.35 1.002 8413.62 2286.62
Between-Person Effects
is_A:pushing_self_cb 0.30 -2.14 2.68 1.001 1750.80 2363.71
is_B:pushing_self_cb 0.93 -1.34 3.37 1.001 2418.81 2811.27
is_A:pushing_partner_cb 0.99 -1.38 3.43 1.001 1999.89 2149.05
is_B:pushing_partner_cb 0.96 -1.46 3.49 1.000 2020.86 2497.22
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb 0.60* 0.03 1.16 1.002 1646.27 2160.85
is_B:barriers_self_cb 0.66 -0.02 1.34 1.001 1654.70 2274.01
Random Effects
sd(is_A) 1.33 0.97 1.82 1.00 1580.22 2402.24
sd(is_B) 1.37 1.03 1.80 1.00 1940.60 2750.37
sd(is_A:pushing_self_cw) 0.13 0.01 0.35 1.00 2190.22 2319.95
sd(is_B:pushing_self_cw) 0.21 0.01 0.52 1.00 1423.00 2184.03
sd(is_A:pushing_partner_cw) 0.28 0.03 0.58 1.00 1546.18 1424.65
sd(is_B:pushing_partner_cw) 0.11 0.00 0.31 1.00 2113.93 2214.85
Additional Parameters
ar[1] 0.24 0.18 0.30 1.00 4694.33 2884.34
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.78 1.72 1.86 1.00 5482.47 2835.78

Compare Models

loo_compare(loo_ind, loo_apim)
##                         elpd_diff se_diff
## pushing_anticipAff_apim   0.0       0.0  
## pushing_anticipAff_ind  -27.0       9.4
bayes_R2(pushing_anticipAff_apim)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2234557 0.008205233 0.2072294 0.2395114
bayes_R2(pushing_anticipAff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2069577 0.008497257 0.1902782 0.2235435

reactance ~ pushing

Indistinguishable

formula <- bf(
  reactance ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day +
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_reactance_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pushing_reactance_ind"
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_reactance_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pushing_reactance_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_reactance_ind)
## Warning: Found 2 observations with a pareto_k > 0.7 in
## model 'pushing_reactance_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
print(loo_ind)
## 
## Computed from 4000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -599.0 24.8
## p_loo        36.9  4.6
## looic      1198.0 49.7
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.7]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     401   99.5%   126     
##    (0.7, 1]   (bad)        2    0.5%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_reactance_ind, ask = FALSE)

summarize_brms(
  pushing_reactance_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.37* 0.15 0.60 1.001 6400.66 3216.87
Within-Person Effects
pushing_self_cw 0.14* 0.02 0.25 1.001 3687.59 2926.71
pushing_partner_cw 0.01 -0.09 0.13 1.001 4397.92 2918.29
day 0.09 -0.32 0.49 1.002 7416.91 3231.80
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw -0.02 -0.09 0.05 1.001 7545.41 3045.25
Between-Person Effects
pushing_self_cb 0.24 -0.33 0.86 1.001 3793.27 3104.16
pushing_partner_cb -0.14 -0.81 0.54 1.001 4135.10 2673.04
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb -0.05 -0.23 0.12 1.000 5844.38 2986.51
Random Effects
sd(Intercept) 0.18 0.01 0.41 1.00 1397.74 2234.57
sd(pushing_self_cw) 0.20 0.07 0.34 1.00 1020.73 1624.63
sd(pushing_partner_cw) 0.10 0.00 0.28 1.00 1736.14 2240.22
Additional Parameters
ar[1] 0.17 0.05 0.29 1.00 4379.90 2715.57
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.02 0.94 1.10 1.00 6018.99 3059.62

APIM

formula <- bf(
  reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 

  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_reactance_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = "models_cache/pushing_reactance_apim"
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_reactance_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with
## `binwidth`.

pp_check(pushing_reactance_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_reactance_apim)
## Warning: Found 5 observations with a pareto_k > 0.7 in
## model 'pushing_reactance_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -601.8 25.4
## p_loo        61.4  7.1
## looic      1203.6 50.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.6]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     398   98.8%   165     
##    (0.7, 1]   (bad)        5    1.2%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_reactance_apim, ask = FALSE)

summarize_brms(
  pushing_reactance_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 0.46* 0.13 0.79 1.002 2786.01 2772.72
is_B 0.31 -0.01 0.64 1.000 4074.33 3123.26
Within-Person Effects
is_A:pushing_self_cw 0.15 0.00 0.31 1.001 2933.18 2937.13
is_B:pushing_self_cw 0.14 -0.02 0.30 1.001 3417.86 3036.79
is_A:pushing_partner_cw 0.06 -0.10 0.23 1.000 2762.82 2567.22
is_B:pushing_partner_cw -0.04 -0.19 0.14 1.001 2832.61 2542.94
is_A:day -0.08 -0.62 0.45 1.000 3794.64 3220.13
is_B:day 0.25 -0.32 0.84 1.001 3897.12 3415.71
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw -0.10 -0.22 0.02 1.001 4940.74 3042.43
is_B:barriers_self_cw 0.05 -0.05 0.15 1.002 5653.14 3029.31
Between-Person Effects
is_A:pushing_self_cb -0.11 -1.19 1.06 1.002 2823.10 2789.67
is_B:pushing_self_cb 0.37 -0.49 1.32 1.001 2738.94 2097.73
is_A:pushing_partner_cb 0.32 -0.89 1.56 1.001 2621.82 2591.60
is_B:pushing_partner_cb -0.19 -1.24 0.96 1.001 2487.53 2450.88
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb -0.04 -0.31 0.24 1.000 2810.45 2892.12
is_B:barriers_self_cb -0.16 -0.50 0.18 1.000 3442.94 3176.45
Random Effects
sd(is_A) 0.40 0.11 0.72 1.00 1224.03 1211.27
sd(is_B) 0.22 0.01 0.59 1.00 1133.35 1771.80
sd(is_A:pushing_self_cw) 0.23 0.03 0.44 1.00 1045.00 861.74
sd(is_B:pushing_self_cw) 0.19 0.01 0.40 1.00 1004.91 1372.28
sd(is_A:pushing_partner_cw) 0.16 0.01 0.45 1.00 1014.73 1886.21
sd(is_B:pushing_partner_cw) 0.15 0.01 0.38 1.00 1451.27 2063.95
Additional Parameters
ar[1] 0.12 -0.02 0.25 1.00 3717.74 3335.13
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.99 0.91 1.07 1.00 3368.06 2991.17

Compare Models

loo_compare(loo_ind, loo_apim)
##                        elpd_diff se_diff
## pushing_reactance_ind   0.0       0.0   
## pushing_reactance_apim -2.8       5.0
bayes_R2(pushing_reactance_apim)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.1937878 0.01849763 0.1570825 0.2293675
bayes_R2(pushing_reactance_ind)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.1609962 0.01983658 0.1206274 0.1984089

Report models for pushing side by side

Indistinguishable Models

pushing_all_ind_models <- report_side_by_side(
  pushing_pa_sub_ind,
  pushing_pa_obj_log_ind,
  pushing_aff_ind,
  pushing_anticipAff_ind,
  pushing_reactance_ind,
  
  model_rows_random = model_rows_random_pushing_ind,
  model_rows_fixed = model_rows_fixed_pushing_ind
) 
## [1] "pushing_pa_sub_ind"
## [1] "pushing_pa_obj_log_ind"
## Warning in summarize_brms(model, short_version = TRUE,
## exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "pushing_aff_ind"
## [1] "pushing_anticipAff_ind"
## [1] "pushing_reactance_ind"
# pretty printing
pushing_all_ind_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_ind()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR pushing_pa_sub_ind 95% CI pushing_pa_sub_ind exp(Est.) pushing_pa_obj_log_ind 95% CI pushing_pa_obj_log_ind b pushing_aff_ind 95% CI pushing_aff_ind b pushing_anticipAff_ind 95% CI pushing_anticipAff_ind b pushing_reactance_ind 95% CI pushing_reactance_ind
Intercept 53.20* [39.52, 70.18] 126.34* [110.67, 143.56] 4.82* [ 4.61, 5.05] 1.85* [ 1.39, 2.30] 0.37* [ 0.15, 0.60]
Within-Person Effects
pushing_self_cw 1.05 [ 0.88, 1.21] 1.03 [ 0.97, 1.09] 0.03 [-0.03, 0.10] 0.00 [-0.14, 0.14] 0.14* [ 0.02, 0.25]
pushing_partner_cw 1.06 [ 0.91, 1.23] 1.03 [ 0.99, 1.07] 0.09* [ 0.02, 0.16] 0.15* [ 0.01, 0.31] 0.01 [-0.09, 0.13]
day 0.66* [ 0.49, 0.88] 0.92 [ 0.81, 1.05] 0.22* [ 0.02, 0.42] -0.25 [-0.71, 0.23] 0.09 [-0.32, 0.49]
weartime_self_cw NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cw 1.10* [ 1.03, 1.17] 1.03* [ 1.01, 1.05] 0.12* [ 0.08, 0.15] 0.18* [ 0.11, 0.26] -0.02 [-0.09, 0.05]
Between-Person Effects
pushing_self_cb 0.95 [ 0.40, 2.07] 0.99 [ 0.70, 1.39] 0.36 [-0.28, 0.96] 0.69 [-0.55, 1.90] 0.24 [-0.33, 0.86]
pushing_partner_cb 0.72 [ 0.30, 1.60] 1.10 [ 0.76, 1.57] 0.23 [-0.43, 0.81] 0.87 [-0.32, 2.10] -0.14 [-0.81, 0.54]
weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cb 1.12 [ 0.96, 1.31] 1.04 [ 0.97, 1.12] 0.27* [ 0.17, 0.38] 0.56* [ 0.30, 0.82] -0.05 [-0.23, 0.12]
Random Effects
sd(Intercept) 0.64 [ 0.43, 0.94] 0.30 [0.22, 0.41] 0.61 [0.47, 0.78] 1.09 [0.82, 1.43] 0.18 [0.01, 0.41]
sd(pushing_self_cw) 0.24 [ 0.03, 0.47] 0.10 [0.02, 0.18] 0.06 [0.00, 0.15] 0.14 [0.01, 0.33] 0.20 [0.07, 0.34]
sd(pushing_partner_cw) 0.21 [ 0.02, 0.45] 0.04 [0.00, 0.10] 0.09 [0.01, 0.18] 0.18 [0.02, 0.37] 0.10 [0.00, 0.28]
Additional Parameters
ar[1] 0.04 [-0.93, 0.94] 0.27 [0.21, 0.34] 0.24 [0.17, 0.29] 0.33 [0.27, 0.38] 0.17 [0.05, 0.29]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.43 [ 0.40, 0.47] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.04 [ 0.00, 0.13] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.53 [0.51, 0.55] 0.86 [0.83, 0.90] 1.85 [1.78, 1.92] 1.02 [0.94, 1.10]
# prepare export to excel
all_models$Pushing_ind <- prepare_for_export(pushing_all_ind_models)

APIMs

pushing_all_apim_models <- report_side_by_side(
  pushing_pa_sub_apim,
  pushing_pa_obj_log_apim,
  pushing_aff_apim,
  pushing_anticipAff_apim,
  pushing_reactance_apim,
  
  model_rows_random = model_rows_random_pushing_apim,
  model_rows_fixed = model_rows_fixed_pushing_apim
) 
## [1] "pushing_pa_sub_apim"
## [1] "pushing_pa_obj_log_apim"
## Warning in summarize_brms(model, short_version = TRUE,
## exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "pushing_aff_apim"
## [1] "pushing_anticipAff_apim"
## [1] "pushing_reactance_apim"
# pretty printing
pushing_all_apim_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_apim()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR pushing_pa_sub_apim 95% CI pushing_pa_sub_apim exp(Est.) pushing_pa_obj_log_apim 95% CI pushing_pa_obj_log_apim b pushing_aff_apim 95% CI pushing_aff_apim b pushing_anticipAff_apim 95% CI pushing_anticipAff_apim b pushing_reactance_apim 95% CI pushing_reactance_apim
is_A 53.36* [36.08, 78.05] 147.66* [126.44, 171.48] 4.92* [ 4.62, 5.21] 1.88* [ 1.31, 2.43] 0.46* [ 0.13, 0.79]
is_B 51.56* [37.38, 70.87] 108.34* [ 93.29, 125.51] 4.68* [ 4.43, 4.92] 1.83* [ 1.24, 2.39] 0.31 [-0.01, 0.64]
Within-Person Effects
is_A:pushing_self_cw 1.02 [ 0.87, 1.19] 1.00 [ 0.93, 1.07] -0.02 [-0.10, 0.07] 0.07 [-0.11, 0.25] 0.15 [ 0.00, 0.31]
is_B:pushing_self_cw 1.09 [ 0.81, 1.38] 1.07 [ 0.99, 1.16] 0.07 [-0.04, 0.17] -0.04 [-0.25, 0.18] 0.14 [-0.02, 0.30]
is_A:pushing_partner_cw 1.11 [ 0.87, 1.38] 1.07 [ 1.00, 1.15] 0.12* [ 0.01, 0.22] 0.16 [-0.06, 0.39] 0.06 [-0.10, 0.23]
is_B:pushing_partner_cw 1.06 [ 0.91, 1.23] 0.99 [ 0.94, 1.05] 0.08 [ 0.00, 0.16] 0.13 [-0.05, 0.30] -0.04 [-0.19, 0.14]
is_A:day 0.58* [ 0.39, 0.88] 0.89 [ 0.75, 1.07] 0.05 [-0.20, 0.30] -0.20 [-0.78, 0.37] -0.08 [-0.62, 0.45]
is_B:day 0.71 [ 0.47, 1.09] 0.98 [ 0.82, 1.18] 0.38* [ 0.14, 0.63] -0.31 [-0.89, 0.27] 0.25 [-0.32, 0.84]
is_A:weartime_self_cw NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cw NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cw 1.06 [ 0.97, 1.17] 1.02 [ 0.98, 1.05] 0.12* [ 0.07, 0.17] 0.10* [ 0.00, 0.21] -0.10 [-0.22, 0.02]
is_B:barriers_self_cw 1.14* [ 1.05, 1.24] 1.04* [ 1.01, 1.07] 0.13* [ 0.08, 0.17] 0.26* [ 0.16, 0.35] 0.05 [-0.05, 0.15]
Between-Person Effects
is_A:pushing_self_cb 0.85 [ 0.18, 3.81] 0.71 [ 0.39, 1.34] 0.83 [-0.54, 2.17] 0.30 [-2.14, 2.68] -0.11 [-1.19, 1.06]
is_B:pushing_self_cb 1.03 [ 0.33, 3.02] 1.10 [ 0.61, 2.01] -0.06 [-1.03, 0.95] 0.93 [-1.34, 3.37] 0.37 [-0.49, 1.32]
is_A:pushing_partner_cb 0.74 [ 0.15, 3.10] 1.43 [ 0.80, 2.56] -0.17 [-1.42, 1.15] 0.99 [-1.38, 3.43] 0.32 [-0.89, 1.56]
is_B:pushing_partner_cb 0.49 [ 0.13, 1.55] 0.95 [ 0.50, 1.79] 0.70 [-0.39, 1.80] 0.96 [-1.46, 3.49] -0.19 [-1.24, 0.96]
is_A:weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cb 1.13 [ 0.84, 1.52] 1.02 [ 0.91, 1.15] 0.21 [-0.08, 0.49] 0.60* [ 0.03, 1.16] -0.04 [-0.31, 0.24]
is_B:barriers_self_cb 1.28 [ 0.99, 1.68] 1.12 [ 0.97, 1.28] 0.28 [ 0.00, 0.55] 0.66 [-0.02, 1.34] -0.16 [-0.50, 0.18]
Random Effects
sd(is_A) 0.82 [ 0.51, 1.26] 0.31 [0.22, 0.43] 0.78 [0.60, 1.01] 1.33 [0.97, 1.82] 0.40 [ 0.11, 0.72]
sd(is_B) 0.55 [ 0.34, 0.82] 0.31 [0.22, 0.43] 0.61 [0.46, 0.81] 1.37 [1.03, 1.80] 0.22 [ 0.01, 0.59]
sd(is_A:pushing_self_cw) 0.11 [ 0.00, 0.34] 0.10 [0.02, 0.18] 0.08 [0.00, 0.21] 0.13 [0.01, 0.35] 0.23 [ 0.03, 0.44]
sd(is_B:pushing_self_cw) 0.34 [ 0.03, 0.80] 0.11 [0.02, 0.21] 0.10 [0.00, 0.25] 0.21 [0.01, 0.52] 0.19 [ 0.01, 0.40]
sd(is_A:pushing_partner_cw) 0.25 [ 0.01, 0.68] 0.08 [0.00, 0.17] 0.12 [0.01, 0.29] 0.28 [0.03, 0.58] 0.16 [ 0.01, 0.45]
sd(is_B:pushing_partner_cw) 0.11 [ 0.00, 0.31] 0.03 [0.00, 0.09] 0.07 [0.00, 0.18] 0.11 [0.00, 0.31] 0.15 [ 0.01, 0.38]
Additional Parameters
ar[1] 0.05 [-0.93, 0.95] 0.22 [0.15, 0.28] 0.15 [0.09, 0.21] 0.24 [0.18, 0.30] 0.12 [-0.02, 0.25]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.44 [ 0.40, 0.48] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.05 [ 0.00, 0.13] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.51 [0.49, 0.54] 0.83 [0.80, 0.87] 1.78 [1.72, 1.86] 0.99 [ 0.91, 1.07]
# prepare export to excel
all_models$Pushing_apim <- prepare_for_export(pushing_all_apim_models)
openxlsx::write.xlsx(all_models, file = "Output/AllModels.xlsx")